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Abstract 



We can lick gravity but sometimes the paperwork is overwhelming - Werner von Braun 

The problem of large scale structure formation in the universe is one of the 
core problems of cosmology today. This thesis discusses some of the issues in- 
volved in explaining how the observed large scale structure in the universe came 
to be. 

This thesis has two distinct parts. The first part (chapters 1-5) discusses the 
issues of structure formation from the view point of standard model of structure 
formation. Chapter 6 discusses alternate cosmologies and structure formation 
scenarios in them. 

In the standard approach we will explore the problem of structure formation 
through gravitational clustering in the universe. The sizes of numbers involved 
renders a microscopic viewpoint cumbersome. We have to adopt a statistical 
procedure which essentially averages over many particles and their individual 
behavior to give us some entities which we shall call "particles" whose behavior 
we can describe more easily. Alternatively we can replace the discrete structure 
of particles by a continuum approximation, thinking of the system of particles 
as a continuous fluid and apply the theoretical models of fluid mechanics to 
describe cosmological systems. In either of the approaches we run into a similar 
problem. When we frame the equations that describe the system, we find that 
the equation has terms which are highly nonlinear, leading to an analytically 

xi 
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intractable system of equations. An analytical solution to the relevant equations 
can be obtained only if one assumes that the system is linear, i.e one drops all 
the terms that are nonlinear in the parameter of interest. But the degree of 
nonlinearity associated with observed structures in the universe like galaxies is of 
the order of 10 3 . Thus we need to deal with the nonlinear terms in the equations 
to establish correspondence with observational results. The nonlinear terms elude 
a simple description, in that a general solution is not yet found. But since we 
have to study the behavior of this equation at highly nonlinear regions we try to 
model the behavior in a variety of ways. 

1. Various approximation schemes such as Zeldovich approximation, frozen 
potential approximation, frozen flow approximation and adhesion approxi- 
mation attempt to take us a little beyond linear theory. 

2. Numerical simulation techniques attempt to integrate the equations either 
in terms of a particle based approach or a field based approach (or a combi- 
nation of the two). N Body simulations form the mainstay of this approach. 

3. An alternate technique treats evolution in time and space as a mapping 
problem and tries to find an appropriate map that takes us from one point 
in the history of the universe to another in a global sense. Scaling laws and 
other ansatzes fall under this category. 

In this thesis we will be concerned with all the the three approaches but with a 
preferential leaning towards the second and third approaches. Herein we explore 
the utility of the third method in various aspects of the study of structure forma- 
tion and the insights it provides into dynamics of gravity in shaping the observed 
structures. 
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The first chapter gives a concise overview of the background. 

Chapter 2 of this thesis uses the third approach to generate nonlinear quan- 
tities from their linear theory counterparts and applies it to the question of "uni- 
versal" profiles in gravitational clustering. We addressed the problem by asking 
the following question: Is it possible to populate the nonlinear universe with 
structures such that the two point correlation function evolves as per linear the- 
ory in all regimes? We find that it is not possible to have strict linear evolution 
but it is possible to find functional forms such that approximate linear evolution 
to any desired order is possible. The earlier investigations had indicated that the 
evolution of density contrast can be separated into three regimes namely linear, 
quasilinear and nonlinear. We have found a functional form which evolves ap- 
proximately linearly in quasilinear and nonlinear ends of the two point correlation 
function. It is also conjectured in this chapter that it should be possible to find 
basis functions which may be based on such approximately invariant forms such 
that the nonlinear density fields can be decoupled to a large extent if expanded 
in terms of these functions. We suggest that the functions that we have derived 
which evolve approximately linearly in all regimes might be a good candidate for 
the role of "units of nonlinear universe" . Another aspect of the analysis tries to 
discover universal aspects of the structures formed via gravitational clustering 
such as density profiles. 

Chapter 3 critically examines the theoretical framework underlying two di- 
mensional N body simulations. If one has to get requisite amount of dynamical 
range in force and mass one requires large grids and large number of particles 
which is not possible given the computing resources available. One attempts to 
get around this problem by trying to do the simulations in two dimensions which 
brings the computational requirements down by a considerable amount. So if it is 
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possible to extract fundamental principles which may be generalized to the case 
of three dimensional gravity from such simulations then they are a good way of 
exploring the fundamental features of gravity. 

There are three ways in which we can define a two dimensional gravitational 
clustering scenario: (i) A set of point particles interacting by 1/r 2 force but with a 
special set of initial conditions such that they are confined to a plane and have no 
velocity components orthogonal to the plane (ii) A set of infinite parallel "needles" 
in which each particle interacts with 1/r 2 force but the "needles" interact with a 
1/r force (iii) A description derived from Einstein's equations in two dimensions. 

Approach (i) is highly contrived and we will not discuss it. The second ap- 
proach based on infinite "needles" suffers from the problem of manifest anisotropy 
because the universe is considered to be expanding in three dimensions while the 
clustering takes place only in two dimensions. In order to explore the third 
approach we developed the formal theory of gravity in D + 1 dimensions and 
considered D = 2 as a special case. The formal analysis of D + 1 dimensional 
gravity led us to the general expressions for scale factor and background density 
in the D + 1 dimensional universe. 

Taking the usual Newtonian limit of the metric and writing down the equa- 
tion describing the growth of density perturbations in the universe via a fluid 
approach made it possible to obtain the D + 1 dimensional analogue of the equa- 
tion describing growth of density contrast. A corresponding formula for spherical 
collapse model is also derived in this work. We then specialize to the case of 
D = 2 and make the following observations. The linearized form of the density 
contrast equations only yield a constant or decaying solution. This is consistent 
with the result that perturbed gravitational potential does not couple to density 
contrast 5. The spherical collapse model solution yields a similar result in the 



1 



XV 



sense that it is not possible to have a gravitational clustering model that grows 
in time. It is possible to obtain clustering by an ad hoc approach by making 
some assumptions but they also lead to inconsistent results in that they give rise 
to singular solutions for the scale factor of the universe. Thus we conclude that 
the infinite "needle" based approach is the only viable way of simulating two 
dimensional gravity. 

In Chapter 4 some of the issues involved in two dimensional gravitational 
simulations are discussed. Nonlinear scaling relations, which have been identified 
in three dimensional simulations, define a mapping from initial linear theory val- 
ues of two point correlation function to the final nonlinear values at a different 
length scale. This allows us to immediately compute the nonlinear parameter at 
a specific length scale knowing the value of the linear one at some other length 
scale. We wish to check the theoretical prediction that similar nonlinear scaling 
relations hold in two dimensions. Another aspect of universal behavior of grav- 
itational clustering that is conjectured is called "stable clustering". This is the 
conjecture that at late times structures have their gravitational infall balanced 
by background expansion leading to a fixed profile. By applying this conjecture 
to the theoretical model for the NSR (nonlinear scaling relations) it is possible 
to derive the theoretical dependence of the NSR in two dimensions if it exists. 
Our conclusions regarding two dimensional gravitational clustering based on the 
simulations are as follows, (i) The prediction is verified and a form of nonlinear 
scaling relation exists for the two point correlation function in two dimensions 
as well. This NSR is independent of the spectrum for the spectra (power laws) 
considered in this study, (ii) In the quasilinear regime the theoretical model 
based on infall onto peaks predicts a dependence that is confirmed by the nu- 
merical experiments, (iii) In the highly nonlinear regime the results diverge from 
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that predicted by the "stable clustering" hypothesis. The "stable clustering" hy- 
pothesis demands that the ratio between infall velocity and expansion velocity 
go to unity but we find that it is driven towards 3/4. Thus we find that stable 
clustering is not a valid hypothesis in two dimensions. 

Chapter 5 addresses the question: What happens at the late stages of clus- 
tering?. The hypothesis of stable clustering asserts that at some point the system 
'virialises'. We examine this process of virialisation and stabilization in more de- 
tail by analyzing what happens to a single spherically symmetric object as it goes 
through the cycle of expansion and collapse. By including a term that describes 
the asymmetries that are generated and enhanced during collapse we examine 
the final state of the system. We begin by writing a modified equation for the 
spherical collapse of a system which includes a term which takes into account the 
asymmetries that are generated during the process of collapse and we derive a 
functional form for this "asymmetry" term by an ansatz that this term depends 
only on density contrast 5. Since we have a relation connecting density contrast 
with the pairwise velocity function h, it is now possible to close the system of 
equations. We require the form of the h function before we can integrate this 
system of equations. By using the fact that the system must reach a constant 
value for h we obtain a functional form for the "asymmetry" term which allows 
us to integrate the equation and analyze the collapse of a single object. This 
leads us to conclude that the system reaches a constant value of 0.65 times the 
maximum radius which is not widely off the value of 0.5 that is usually stipulated. 
Thus we demonstrate that the growth of asymmetries can be used to stabilize 
the collapse. 

In Chapter 6 the question of approaching the problem of structure formation 
from other alternate cosmological scenarios is discussed. To address this ques- 
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tion in the context of Quasi Steady state cosmology is the theme of the work 
discussed here. An algorithm was devised which mimicked the basic physical 
content of QSSC (Quasi Steady State Cosmology) model by the following geo- 
metrical method. We generate N random points and require that in the next cycle 
the volume becomes eight times the initial volume (due to the expansion) and the 
particles generate new mass particles in their vicinity. Then we select the central 
volume equal to the initial volume and repeat this scale and shrink process. To 
apply this algorithm in the context of QSSC model we use the fact that only a 
fraction of the particles given by / = 3Q/P where Q and P are parameters of 
the QSSC model create new masses in the next generation. Consequently the 
scaling used is exp(/). We have succeeded in showing that this model reproduces 
the observed two point correlation function with a slope of —1.8. Visual analysis 
shows clear indication of the clustering growing and the initial smooth distribu- 
tion separating into voids and clumps. This approach is characterized by the fact 
that clustering in this model takes place without the help of gravity, i. e gravity 
plays no role in inducing clustering. 

The final chapter of the thesis presents the conclusions and integrated dis- 
cussion based on the previous chapters. 
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Chapter 1 
Introduction 



Duct tape is like the force. It has a light side, a dark side 
and holds the universe together. - Carl Zwanzig 

Observations of the universe in multiple wavelengths have revealed large ag- 
gregates of matter on all length scales. Figure 1.1 shows the distribution of 
galaxies from Las Campanas Redshift survey which clearly shows clustering in 
the distribution of galaxies. The formation and existence of large scale structures 
such as galaxy clusters and superclusters is an important problem in cosmology. 
In spite of the inhomogeneities on these scales it is also observed that on very 
large length scales (> 300 Mpc) the universe appears homogeneous. This permits 
the universe to be modelled as being made up of statistically similar volumes 
of linear dimensions of about 300Mpc or more. Standard cosmological models 
based on General Theory of Relativity are derived from this assumption of sta- 
tistical homogeneity and isotropy of the universe. The problem of large scale 
structure formation in these models thus involves explaining how the observed 
inhomogeneities and anisotropics were initially generated in this uniform back- 
ground and how they grew into the observed structures. Once a mechanism for 
generating the initial perturbations of the smooth distribution of matter is pos- 
tulated the transition from uniformity on large scales, to the highly non uniform 
structures at small scales, can be explained by gravitational clustering. Obser- 
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Figure 1.1: The Las Campanas Redshift survey (http://www.lcrs.com) 

vations by space experiments (C0BE,WMAP etc.) have revealed the presence 
of small inhomogeneities in the uniform density fields of matter at early times. 
It is assumed that these seeds evolved via gravitational clustering into the ob- 
served large scale structures and experienced a multitude of non gravitational 
phenomena at small scales (such as star formation) to evolve into the present day 
universe of clusters and galaxies. 

This chapter introduces the basic framework of this model in brief. It is orga- 
nized as follows. The first section 1.1 summarizes the main features of the smooth 
background cosmological models. The following section 1.2 discusses the equa- 
tions governing the dynamics of matter in this expanding background cosmology 
in discrete and continuum approximations. Section 1.3 discusses the solutions of 



this nonlinear system of equations in linearised regime. Nonlinear approximation 
schemes such as Zeldovich approximation and spherical top hat approximation 
are discussed in section 1.4. Section 1.5 deals with N body computer simulations 
as a way of modeling and studying some features of nonlinear structure forma- 
tion. Various statistical measures used to quantify and analyze the formation of 
structures are discussed in section 1.6. Section 1.7 deals with an alternate ap- 
proach to time evolution of some statistical quantities, thus obviating the need 
for expensive and time consuming full scale N body simulations. 

1.1 Background Cosmology 

A model appropriate for a universe dominated by gravity on large scales, is best 
formulated in the context of General theory of Relativity. This model, which 
will be referred to as the 'background cosmology', is described by the metric 
of the spacetime.The evolution of the metric is governed by the energy density 
distribution through Einstein's equation given by 

G\ = R k -UlR = MGT k (1.1) 

where G\ is the Einstein tensor, G the gravitational constant and T % k the energy 
momentum tensor. The 'cosmological principle' which assumes a statistically 
homogeneous and isotropic universe is further used to constrain the form of the 
metric of the universe to ([36], [38], [45], [46] ) 

dr 2 



ds 2 = dt 2 -a 2 it) 



+ r 2 (d9 2 + sin 2 6d(f) 2 ) 



(1.2) 



1 — kr 2 

(we use the convention c = 1 consistently throughout unless explicitly indicated). 
This metric is completely characterized by a time dependent function scale factor 
a(t) and curvature k. The constant k determines the geometry of the universe - 



(the curvature of the spatial hypersurfaces of the Friedmann universe) — taking 
the values k > in a closed universe, k = in a flat universe and k < in 
an open universe. The 'fundamental observers' to whom the universe appears 
isotropic and homogeneous remain at constant (r, 9, 0) with their physical sepa- 
ration increasing in proportion to scale factor a(t) which determines the overall 
scale of the spatial metric. An observational consequence of expansion is that 
the light emitted by a source at time t is observed now, at t — with a cos- 
mological 'redshift' z = a /a(t) — 1 where a = a(t = 0). Another quantity 
of primary significance in observational cosomology is the 'Hubble Constant' H 
defined by H = a/ a, which measures the rate at which the universe is expanding 
at a given point in time. The value of 'Hubble constant' today is denoted by 
H = lOOhkms' 1 Mpc" 1 with 0.52 <h< 0.62 [55]. 

The 'cosmological principle' also constrains the source term in Einstein's 
equations, namely the energy momentum distribution, so as to make it con- 
sistent with the assumed homogeneity of the metric, leading to a diagonal form 
Tp = dia[p(t), —p(t), —p(t), —p(t)]. This form is motivated by the assumption of 
the source to be an ideal fluid with pressure p and density p, whose stress tensor 
has the above form in the rest frame of the fluid. Under these assumptions, an 
extra equation of state connecting p and p allows the source to be completely 
determined. 

When this source term is plugged into equation 1.1, the Friedmann equations 

a 2 + k 8ttG 



a 2 3 
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1 5— = -8nGp (1.4) 

a a z 

for the scale factor a(t) which determines the metric given the curvature k, are 
obtained. 



From equation 1.3 we can define the value of curvature k by 
k 8ttG 



Oq 3 



-po - Hi = H 2 (Q - 1) (1.5) 



where ao is the scale factor at present and Q = po/p c where po is the present 
value of density and p c = 3Hq/8itG is called the critical density. Equation 1.5 
shows that the universe will be open, closed or flat depending on the values of fl 
being less than, more than or equal to one respectively. 

We assume that the total p of the universe may be separated into contribu- 
tions from the various constituents of the the cosmic fluid. Other than radiation 
and baryonic matter there is compelling evidence for the presence of a dominant 
component that does not interact with light and is consequently not directly ob- 
servable, called the 'dark matter'. There is also good observational evidence for 
the presence of a vaccum energy density p v [3]. For each component of energy 
density p with an equation of state given by p = p(p) , the density varies with the 
scale factor according to energy conservation 

d(pa 3 ) = -pd(a 3 ) (1.6) 

For a generic equation of state given by p = wp, equation 1.6 gives p oc 
a -3(i+i«)_ ^he e q Ua tions of state for radiation component of the fluid is p = |p 
and for the non relativistic pressureless dust, p = [36]. Defining Q x = p x /p c 
where the subscript x is r for radiation, nr for nonrelativistic matter which has 
contributions from both baryonic matter p B and 'dark matter' Pdm, v for the 
vaccum energy density, we can write the complete time dependence of the scale 
factor a(t) as 



Ptot{a) — Pcrit 



« 2 TT2 

a 2 
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The radiation density varies as a~ 4 whereas the matter density goes as a -3 . 
This has the interesting consequence that although radiation density is much 
smaller than matter density today, in the past the radiation field dominated. The 
point in time where the densities were equal is defined as the 'matter-radiation 
equality'. The dynamics of the universe from this epoch is dominated by matter, 
and we speak of 'radiation dominated epoch' as opposed to 'matter dominated 
epoch'. The redshift of equality z eq and the value of Hubble constant at equality 
H eq are defined by, (neglecting curvature and radiation components) 



(X e q Q r 



[l + z eq ) (1.9) 



H 2 eq = 2H 2 Q nr (l + z eq f (1.10) 

Since the curvature term and the vaccum energy terms do not contribute signifi- 
cantly in early times [36], analytical solutions to 1.7 valid for t ^> t eq and t < t eq 
may be obtained, where t eq is the time of equality. 

/ o \ 2 /3 

" - ' — I <H eq tf' z (fart^t eq ) (1.11) 
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Figure 1.2 shows the complete evolution for the scale factor a(t) (obtained 
by integrating Eq.(1.8)with Q v and (fl — 1) terms neglected) with the solutions 
in the matter dominated and radiation dominated regimes superposed. 

In the standard model,the initially hot universe whose evolution is dominated 
by radiation energy transitions to a matter dominated structure at z ~ 10 4 . In 
this hot universe, electrons and photons couple to each other via scattering pro- 
cesses such as Thompson scattering. At z ~ 1100 the temperature cools enough 
for neutral hydrogen to form and the photons and matter decouple. The photons 
then travel freely, cooling adiabatically, until they are observed at the present time 
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Figure 1.2: The scale factor (solid line) for a f2 = 1 universe with the solution for 
radiation dominated epoch (dashed line) and matter dominated epoch (dot dashed 
line) superposed on it. 

as the 2.73K cosmic microwave background (CMB) radiation. The observations 
of CMB ( [1],[2]) show that the universe at the time of formation of neutral hydro- 
gen (last scattering surface) was extremely uniform with fractional fluctuations 
in energy density ( and consequently in gravitational potential) of approximately 
10~ 5 . These small fluctuations, are the seeds of the large scale structures seen 
today. They are amplified through gravitational interaction, eventually leading 
to the formation of galaxies and stars. 
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1.2 Dynamics of matter 

To analyse the evolution of matter in this expanding background a dynamical 
model which describes how matter aggregates under the influence of force of 
gravity in an expanding universe is required. 

If the system has to be described in it's full generality one must use a general 
relativistic description. But one can use a 'Newtonian' approximation valid for 
certain regimes by finding the effective 'Newtonian limit' of the Friedmann metric 
[38] . Applying the transformation to new variables R and T defined by 

R = ra(t) (1.13) 

T = t - t + -adr 2 + C(r 4 ) (1.14) 
2 

on eq 1.2 we get 

ds 2 » {l-'-R 2 )dT 2 -{l + ^-R 2 + ^-R 2 )dR 2 - 

a a 2 a 2 

R 2 {d6 2 + sin 2 6d<f) 2 ) (1.15) 

Near R = one can assume a locally inertial coordinate system if one restricts to 
quadratic order in R/dn where dn = cH^ 1 is the 'Hubble Radius'. In the weak 
field limit, we know that goo = (1 + 20a?) where <pN is the Newtonian potential. 
Hence the equivalent Newtonian potential of the FRW metric is 

<i>FRw(R,t) = -\~R 2 (1.16) 

To proceed further one has to investigate how particles move in a self consistent 
manner in universe described by a perturbed Friedmann metric, with the per- 
turbations introduced by the potential of the particles. Since in the limit of 
weak gravitythe perturbed metric can be written as a linear superposition of the 
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effective gravitational potential of the unperturbed metric and the perturbation 
due to matter, the effective Newtonian potential of the perturbed metric is 

<t>N = <f>FRW + <t> I 1 - 17 ) 

In the transformed coordinate system we can write the equation of motion 
for a particle in the Newtonian limit as 

^ = ^(«x,) = -V^ (1.18) 

where Xj is the coordinate of the particle i. This equations when expanded, 
and the contribution from the effective Newtonian potential of the FRW metric 
cancelled, leads to 

d 2 1 
a— x, + 2ax, = — V x <j) (1.19) 
dl z a 

where is the perturbed potential which satisfies Poisson's equation 

Vl(j) = ArfGa 2 p bm 5 (1.20) 

Pbm is the smoothed density of the background and 5 = (p/ pb m ) — 1 is the density 
contrast which defines the perturbations on the smooth background. 

To the same order accuracy, we can replace d 2 /dT 2 by d 2 /dt 2 in the equation 
of motion leading to 

x, + 2-Xi = -^V x (f) (1.21) 
a a 2 

as the equation of motion governing the trajectory of a particle moving in an 
expanding universe in a potential created by the perturbation 5. 

The density field p(x, t) of a set of point particles may be defined by the 
following expression 

in 

p(x) = ^5 Dirac (x- Xi (t)) (1.22) 

a \ l ) i 

(1.23) 
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Taking derivatives, Fourier transforming to get the density modes in Fourier space 
5 k (t) and using the equation of motion for x we get [38] 

k + 2-5 k = A7iGp bm 5 k + A k -B k (1.24) 



where the A k and B k terms are given by 



k' 2 



kVO.k 

B k = ^^(k-Xjfexp^k-x^)] 



k • k' k • (k - k') 

+ 



Ik-k 



'12 



(1.25) 
(1.26) 



The two terms A k and B k are the nonlinear mode coupling terms which cause 
the evolution to have a nonlinear structure. The summation in term B k is over 
all the particles in the system. Neglecting A k and B k when 5 is small leads to a 
linear equation that has a solution that grows as well as a solution that decays 
in time. 

The equations that govern the mechanics of the fluid in the continuum are the 
continuity equation, the Euler equation and the Poisson equation which governs 
the evolution of the potential. 



^ + V r -( Pm U) = (1.27) 
ot 

^ + (U-V)U = -V0 tot = -V0 FW -V0 (1.28) 
ot 

V 2 r <j ) = 47rG(p m - Pbm ) = AnG Pbm 5 (1.29) 

Using a different time parameter b(t) defined as growing mode solution of the 
linearised equation for density contrast and the corresponding peculiar velocity 
u = v/(ab) where v = ax is the original peculiar velocity, results in a second 
order equation which governs the growth of density contrast 5 [38] 

d 2 5 3Ad5 3A C/ rx 4 1 ( al5\ 2 . PW 9 ^ 
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in the fluid approximation, a and Q represent the shear and rotation of the 
peculiar velocity field of the fluid given by o 2 = <7 o6 <7 a 6 and Q 2 = Q a Q a . and 
A = (p bm /pcrit)(ab/ab) 2 . 

Solutions to these formal mathematical equations are rather difficult except 
in special cases, due to the presence of the highly nonlinear terms. Understanding 
the nonlinear terms is the key to understanding structure formation, since many 
of the observed structures are highly nonlinear entities with the parameters of 
interest such as density contrast of the order of a thousand or more. 

1.3 Linear perturbation theory 

Solutions to the equation governing the growth of 5 may be obtained in the 
framework of linear perturbation theory [36] . In this approach the nonlinear terms 
A k and B k are assumed to be zero due to the small values of 5 and velocities. 
The equation that governs the growth of 5^ in the linear regime is 

5 k + 2-5k = 4nGp bm 8 k (1.31) 
a 

This equation as discussed before yields solutions which describe growth as 
well as decay of density contrast. The growing mode solution shows that 5 in 
matter dominated phase evolves proportional to scale factor a in an Q = 1 uni- 
verse. The solutions for time dependence of potential and velocity field indicate 
that the gravitational potential remains constant in time and the velocity field 
in linear regime is proportional to the gradient of potential given by 

where / ps f2°- 6 . 

Since the regimes of interest are highly nonlinear, we require insight into 
the behaviour of the nonlinear terms. One of the important physical effects 
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that can be readily mapped into the nonlinear terms is obvious from equation 
(1.25). This equation clearly shows that mathematically speaking the nonlinear 
terms couple various Fourier modes of density contrast field in that it causes 
evolution of one mode to be dependent on the other. This gives an equivalent 
way of identifying the linear regime as the period of evolution of a mode, where it 
evolves independently of the others. The slow growth of nonlinearity causes this 
independence to be broken until the strong coupling at highly nonlinear scales 
links together many modes of the system. 

1.4 Nonlinear approximations 
1.4.1 Zeldovich approximation 

The approximation schemes attempt to either extend the linear theory results 
such as constancy of potential, into nonlinear regimes or model the universe by 
analytically tractable structures such as spheres. The most fruitful of the ap- 
proximation schemes of the first kind is the Zeldovich approximation [58]. This 
is a simplistic yet effective approximation which assumes that the particles move 
in inertial trajectories, with their initial velocities given by the initial potential. 
These trajectories cross, leading to caustic surfaces where high density aggregates 
of matter are produced. The shapes of the caustic surfaces of the initial poten- 
tial/density field determine the shapes of the large scale structures, which are 
formed in this approach. Thus this approximation scheme is a simple one time 
map from Lagrangian space q to the Eulerian space x, given by 





where b(t) is the growing mode solution from linear perturbation theory and f is 
the initial velocity field oc V0. The density in Eulerian space 

^ = [(1 - bit) a) (1 - bit) 0) (1 - b(t) 7)]- 1 (1.34) 
Po 

where (—a, —(3, —7) are the eigenvalues of <9 2 f / dq 2 . This expression clearly shows 
that the collapse will take place along the axis defined by the largest negative 
eigenvalue. The distribution function of the eigenvalues for a gaussian random 
field indicates that collapse along one dimension in that has a much higher prob- 
ability of occurring since a 3> f3, 7. Thus the first structures that form will be 
sheet like structures that are likely to be highly warped. 

The following figures (Fig. 1.3 - 1.6) show a slice of a full N body simulation 
and a corresponding scenario generated by the Zeldovich approximation. It is 
visually evident that Zeldovich approximation is fairly accurate. The point of 
breakdown of this approximation comes when the particles after having reached 
the caustics start moving away from it, i.e the shells start crossing and passing 
through each other leading to a thickening and dissolution of the structures that 
form (Fig 1.6). This approximation when supplemented by an artificial viscosity 
term (the 'adhesion approximation') is more faithful to the N Body simulation 
even at late times. There are many variants on the basic Zeldovich approximation 
such as the truncated Zeldovich approximation which are all consistent with N 
body simulations to various degrees. Some of the other approximation schemes 
driven by physical considerations are frozen potential approximation, frozen flow 
approximation and so on which extend the linear theory results related to density, 
velocity and potential field growth in various ways. Nonlinear growth of a single 
object is also of interest because the universe at late stages can be modelled as 
a set of highly collapsed objects whose overall configuration in the universe is 
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Figure 1.3: N-body simulation at a=0.25 (n=-l. Power law spectrum). 

determined by the large scales in the power spectrum of perturbations. Spherical 
collapse is an approach in this direction. This model has been successfully used 
to derive evolution of mass spectrum, non linear scaling relations and so on. 

1.4.2 Spherical model 

A single virialised structure such as a galaxy or a cluster in the universe is mod- 
eled as a spherically symmetric density perturbation in this approximation. The 
assumption of spherical symmetry allows the system to be analytically solved 
[4] and the whole history of evolution from initial conditions till the collapse 
and stabilization to be dealt with analytically. The system is modeled as con- 
sisting of concentric spherical shells of radius Ri(t). It can be seen that this 
overdense/underdense sphere can be modeled as an embedded closed/open uni- 
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Figure 1.4: Zeldovich approximation at a=0.25 (n=-l. Power law spectrum). The 
correspondence between structures in the full Nbody simulation is evident. 

verse and the equation of motion for the shell can be identified with the equation 
of growth of scale factor. This analogy permits the solution for time evolution of 
a shell to be written in a parametric form as 



Ri = A(l - cos 9) 
t = B(6- sin6) 



(1.35) 
(1.36) 



where A 3 = GMB 2 (M is the total mass inside Ri). Defining the density contrast 
in terms of the average density within the sphere S, its time dependence may be 
approximated as 

3 /6t\ 2/3 

5~-(-) (1.37) 
20 \BJ v ; 

Comparing with the linear theory results in an Q = 1 universe it can be seen 
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Figure 1.5: N-body simulation at a=0.5 (n=-l. Power law spectrum). 

that the sphere will breakaway from the general expansion and reach a maximum 
radius at 6 = n. At this point the linear theory result for density contrast is 
5un = 1.06 where as the nonlinear density contrast is 5 sp h = 5.6. The final fate of 
the sphere will be an uninterrupted collapse towards a single point at 6 = 2n. The 
linear density contrast at this time is 5u n = 1.68 where as the actual nonlinear 
density contrast essentially shoots up towards infinity (Figure 1.7). 

To prevent the model from being driven towards a singularity one has to 
invoke a physical mechanism referred to as 'virialisation'. This involves an ad 
hoc principle of introducing a stability criterion such that when the potential 
energy is equal to twice the kinetic energy the collapse will stop and the system 
will reach a final density contrast of about 180. 
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Figure 1.6: Zeldovich approximation at a=0.5 (n=-l. Power law spectrum). The 
thickening and dissolving of the structures formed is clear leading to breakdown of 
Zeldovich approximation. 

This approximation has proved to very effective, despite its highly restrictive 
assumptions, in quantitatively describing measures such as mass functions and in 
theoretically modeling nonlinear scaling relations. However spherical symmetry is 
a restrictive assumption and a model which is more representative of the observed 
structures is to be sought. 

1.5 N Body Simulations 

An effective model of the phenomena so that the process can be better understood 
must have the following three components (1) a representation of the basic entities 
that are being modeled (2) a set of initial conditions (3) a dynamical system of 




Linear density contrast 



Figure 1.7: The nonlinear density contrast versus the linear density contrast for the 
spherical collapse model 

equations to evolve the system. The various models differ in the ways in which 
these three requisites are met. Truly speaking, they also differ in the domain of 
the problem that they address but it can be traced down to the way in which 
they implement the above conditions. 

(f) The basic entities must be chosen such that the high degree of detail af- 
forded by computer simulations is not detrimental to understanding these struc- 
tures better. A complete model of the universe taking into account all the funda- 
mental interactions between the particles is therefore ruled out. A more feasible 
and productive approach is to treat large mass aggregates as the "particles" of 
the system and describe the universe in terms of the dynamics of these massive 
"particles" which are governed by gravity. This approach forms the basis of N 



body simulations. An alternate approach involves averaging over the local veloc- 
ities of the particles in a small volume and treating the volume itself as a fluid 
element and describe the system in terms of dynamics of a fluid. Analytic ap- 
proaches usually favor this picture. Observationally all the matter in the universe 
consists of that which is inferred from the various signals that we receive via the 
electromagnetic spectrum. Some observations such as rotation curves of galaxies 
as well as velocities of galaxies in clusters indicate that there is a large, dominant 
component of matter in the universe, that has no interaction with light. This 
matter is referred to as 'dark matter'. Many observational problems such as the 
observed flat rotation curves of galaxies and so on can be explained by a judi- 
cious invocation of this invisible "matter" whose constituents and interactions 
other than gravitational, is yet to be discovered. Since it is indicated that more 
than 90% of the observed universe consists of dark matter,the models incorporate 
the dynamics of the dark matter prior to that of visible matter. Thus the com- 
puter simulation model has three core components. The Friedmann metric, the 
dark matter and the visible matter. The dark matter forms the major component 
as far as gravitation is concerned. 

The various forms of dark matter that have been used are (1) CDM [Cold 
dark matter] (2) HDM [Hot dark matter] and a mixture of the two. The hot dark 
matter and the cold dark matter differ in the way power is distributed in various 
modes. The most notable physical difference is due to 'free streaming' which 
causes small scale power to be suppressed in the case of hot dark matter. This 
qualitatively and quantitatively changes the way hierarchies of structures form in 
the two models. The next set of entities which are baryons (or the components 
of visible matter) are dealt with separately since their gravitational contribution 
is much smaller than that of dark matter but they have interactions which are 



non gravitational which govern their dynamics. Thus the basic entities and their 
interactions are as follows. The dark matter gravitationally clusters in an universe 
with an expanding background. The baryonic particles collect in the deep dark 
matter potentials with their structure and state dominated by additional non 
gravitational processes such as heating, cooling and so on. 

(2) The initial conditions: The universe that is observed is considered to be 
a single realization of an underlying gaussian random field which is predicted 
by various theories of early universe and has been observationally supported by 
the COBE data. The gaussian random field is wholly described by the power 
spectrum defined as follows. 

P(k,k') = |<5 k | 2 (1.38) 

where <5k) is kth Fourier mode of random field 5. 

The quantities that are used to describe the structure of the universe are also 
statistical in nature. This implies that exact structure of the universe as observed 
today cannot be modeled by these approaches, but rather allows us to describe 
a statistical model of the universe. 

(3) The dynamics; The dynamical equations that describe the growth of the 
system when it is treated as a set of interacting particles consist of the equation of 
motion for the trajectory of a single particle in an expanding background and an 
equation which allows us to evaluate the force on a particle due to all other parti- 
cles. N Body simulations essentially define a number of particles and generate a 
phase space configuration for them based on the chosen initial conditions which 
are usually specified by the initial cosmological model, ie the parameters that 
describe the background universe such as density of matter, density of vacuum 
energy, Hubble's constant etc and the power spectrum of the initial perturbation 



field. The initial conditions also specify the epoch at which these conditions are 
given. This epoch is usually chosen to be the epoch of recombination at a redshift 
of approximately 1000. 

The initial conditions can be evolved analytically to a point where the non- 
linear effects start to increase in importance and from there the system is allowed 
to evolve according to the full system of equations until the desired point in time. 
This approach is detailed in the following review [9] . 

N Body approach is complete but not feasible to a high degree of accuracy 
owing to the computational limitations that exist. The limitations are usually 
caused due to resource (memory and disk storage) and computational (CPU 
power, time) requirements being inadequately met by existing technology. The 
large number of theoretical models and the continually increasing accuracy of 
observational data require that the system be modeled more and more accurately. 
In addition the outputs of simulations which consist of the phase space of all the 
particles in the system are too detailed and a large amount of further computation 
has to be performed before statistical measures are derived from them which may 
be used to effect the comparisons. 

These limitations led to attempts to model the system under study by tech- 
niques which either bypass the simulation strategy (analytic and semi analytic 
techniques) or use some physically motivated ansatz to arrive at a final config- 
uration from an initial one, essentially by bypassing some of the steps in the 
simulation. The second category of approaches include the various numerical 
approximation schemes such as Zeldovich approximation, frozen potential ap- 
proximation and so on. The other approach leads to perturbative techniques 
for exploring the quasi and nonlinear regimes as well as some nonperturbative 
techniques that attempt to look at the problem of evolution in time as a one 



shot time dependent mapping which gives some of the statistical measures com- 
putable from the final configuration in terms of the measures computed over the 
initial configuration at one go. Thus the three approaches, the complete N body 
approach, the approximation schemes and the mapping approach, all attempt to 
explain the formation of structure in a complimentary sense. 

The algorithms for N body simulations have been discussed in the references. 
The different categories of simulations are essentially based on the force compu- 
tation method. This thesis will be primarily using a PM (particle mesh) code 
developed by Bagla and Padmanabhan [9]. The N body experiments and some 
of their results are discussed in this section. There are two kinds of PM N body 
simulations discussed here: the two dimensional simulations and the three di- 
mensional simulations. The 2D simulations allow us to explore a much higher 
dynamic range in the given computational resources and consequently many ex- 
periments are conducted in two dimensions. The 3D simulations correspond to 
full scale models for the universe at hand. The essential nature of gravitational 
clustering in terms of how the "particles" of the system respond is easily visu- 
alized in two dimensions. The following results reveal that initially the particles 
are distributed in a random manner according to chosen power spectrum. Then 
the clustering proceeds forming filamentary structures whose intersections serve 
as the points to which matter flows along the filaments to form clusters. In three 
dimensions this phenomenon is essentially repeated with sheets, filaments and 
clusters. 

Figures(1.8 - 1.10) which show the results visually for a set of 2D simulations 
clearly reveal the initial peaks and the formation of long filamentary structures 
and the later collapse into the peaks themselves. To the first approximation 
one can clearly understand the fact that the initial peaks in the system serve 




100 200 300 400 

Figure 1.8: The flow field at a=0.5 



as the nucleation centres for structure formation and these nucleation centres 
themselves do move due to the power existing at length scales corresponding 
to their mean separation. The question of first structures that form, whether 
they are filamentary or sheet like in 3D and it's analogue in 2D has partial 
answer in terms of the Zeldovich approximation. This approximation scheme 
which essentially places the particles in inertial motion with their initial velocities 
given by the initial conditions indicates that the initial structures are sheet like 
(string like in 2D) which is roughly concordant with what the simulations indicate. 
Observationally too we see the presence of very large structures (100 Mpc) which 
may be sheet like or filamentary. A simple analysis based on the Zeldovich 
approximation which is a valid approximation until shell crossing takes place 
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Figure 1.9: The flow field at a=1.5. It can be seen that the velocity vectors flow 
along the caustics (filaments) towards the cluster centers 
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gives the statistical answer, that for a gaussian random field it is more probable 
that sheets would form rather than filaments in 3D and correspondingly in 2D. 
The pictures below clearly reveal the initial formation of filamentary structures 
(the arrows indicate the direction of velocity) and the damping of the velocity 
transverse to the filament. This causes the matter to flow along the filaments 
into the junctions where filaments meet which are essentially where the roughly 
spherically symmetric clusters form. This picture is expected to hold good in 
three dimensions also. 
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Figure 1.10: The flow field at very late times. The flow is almost completely concen- 
trated at the centers. The particles participating in these flows transfer power from 
large scales to small scales 

1.6 Statistical indicators 

In developing the framework discussed in the earlier sections, we have made 
some simplifying assumptions regarding homogeneity, Newtonian limit of the 
metric and so on. Since these approximations are valid only at different length 
scales (homogeneity on large length scales and Newtonian approximation on small 
length scales) it is clear that our models for the universe as a whole is built from 
smaller models at different length scales. To extend the results derived from 
the smaller models to the universe as a whole as well as correlate them with 
observations, we need to use statistical approaches. Thus the models allow us to 
predict various statistical quantities which may be compared to observations. 



27 



The statistical quantities that are used to analyse and compare the N body 
data with other simulations (as well as with observations) are computed from 
the positions and velocities. The statistical quantities are usually measures of 
deviation from the assumed smoothness and distributions of such quantities such 
two point correlation function, power spectrum and higher moments of the dis- 
tribution functions. Another approach followed when models are being compared 
to observational data, is to constrain the parameters of the model by using statis- 
tical goodness of fit and other tests of significance. In this section we shall briefly 
discuss some of the quantities relevant to later discussions. 

The underlying gaussian field is completely characterised by the power spec- 
trum P(k) defined by 



One can further define other related quantities such as two point correlation 
function in terms of the power spectrum 



The assumed statistical isotropy of the universe implies that the two point cor- 
relation function and power spectrum are functions of magnitude of r and k 
respectively and is independent of direction. Some of the other measures of 
departures from homogeniety are the variance of the density field after an appro- 
priate smoothing using a window function such as 



P(k) = \5 k \ 2 = P(k) 




(1.39) 




(1.40) 
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[sin kR — kR cos kR] 



(1.41) 



k 3 R 3 




(1.42) 
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giving 




(1.43) 



(1.44) 



Figure 1.11 shows the quantities discussed above computed for a standard COBE 
normalized CDM power spectrum given by 



where the parameters A, B,C, D,n are equal to (24/h) 4 Mpc 4 , 1.77 /(Qh 2 )Mpc, 
^tth 2 )- 1 - 5 Mpc 3 / 2 , (tth 2 )- 2 Mpc 2 and 1 respectively (h = 0.65)[36]. 

A more complex quantity is the Press Schecter mass function (Fig. 1.12) which 
gives the number of collapsed objects of mass M as a function of redshift 



As has been indicated, all these different measures which are low order mo- 
ments of the underlying distributions are not independent of each other. Higher 
order moments contain more detailed information which is required to discrim- 
inate between models for instance when the lower order correlations are unable 
to do so. There exist statistical measures such as minimum spanning tree, per- 
colation statistics and so on which attempt to deal with the whole hierarchy of 
moments in a holistic manner. They can be used to characterise complementary 
aspects in the statistical description of large scale structure. Various other quan- 
tities are also used to quantify morphologies and topological properties of large 
scale structure such as minkowski functionals which are related to the moments 
of the distributions. 
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Figure 1.11: Different statistical indicators such as the two point correlation function 
and A 2 (k) = k 3 P(k)/(27i 2 ).(A(k) - dashed line, £ - solid line, a sph - dotted line, 
&gauss ~ dash dot line for a standard CDM power spectrum with k = 1/7 

1.7 Nonlinear scaling relations 

The N body approaches form an expensive way of arriving at the statistical 
properties of the late universe given the statistical quantities describing the early 
universe. An alternate approach which is a one time mapping approach to nonlin- 
ear evolution is referred to as 'Nonlinear scaling relations'. It was demonstrated 
empirically from simulations by Hamilton et al [19] that the averaged two point 
correlation function in the nonlinear regime may be obtained from the initial 
average two point correlation function via a non local map given by 



Inl{x) = f NL {L(l)) 



(1.47) 
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Figure 1.12: The mass function plotted as the log of fraction of mass in struc- 
tures more massive than M (ln(f2)) as a function of a(M) for different redshifts 
z=0, 1,2,3,4,5 from top to bottom. 

where x and I are related by / 3 = (1 + ^nl(x))x 3 . The functional form for f NL 
proposed by Hamilton et al [19] is as follows 



It is possible to give a theoretical explanation for the existence of such a mapping 
[34] through an argument based on the pair conservation equation which leads 
to the result that the universal mapping function for the two point correlation 
function may be obtained from the behaviour of h — —v/(ar), ratio of average 
pair velocity to Hubble velocity ar under an assumed closure condition. This 
allows us to express the mapping function as an integral over the h function. 
A theoretical modeling of the mapping based on infall on to high peaks in the 



z + 0.358 z s + 0.0236 z' 



(1.48) 



1 + 0.0134 z 3 + 0.0020 z 9 / 2 
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intermediate regimes [37] allows the empirical relation to be well approximated by 
a discontinous powerlaw approximation with three different values for the power 
law indices in the three regimes, linear, quasilinear and nonlinear. The power law 
index in the nonlinear regime is related to the value of h in the highly nonlinear 
end which might not be unity as is implied by the stable clustering hypothesis. 
This approach permits the notion of a more generalised stable clustering to be 
introduced with a corresponding slope for the nonlinear end of the mapping. The 
power law representation in three dimensions is 



Jnl{z) 



' 1 (z«l) 
z 3 (1 < z < 200) 
z 3 ' 2 (200 < z) 



1.49) 



Figure 1.13 shows the mapping function in three dimensions and power law 
approximations to the same. Figure 1.14 shows the evolution of average two point 
correlation function for a CDM power spectrum at a late epoch and the linearly 
evolved two point correlation function at the same time. The nonlinear evolution 
of the two point correlation function may be seen at small scales. 
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Figure 1.13: The nonlinear scaling relation between nonlinear and linearly scaled two 
point correlation function. 
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Figure 1.14: The linearly scaled two point correlation function (solid line) and the 
nonlinear two point correlation function of CDM power spectrum at a late epoch 
(a=l) 

In brief the approach taken towards studying structure formation in the uni- 
verse is a) a model for the underlying smooth background universe b) a model for 
the seed perturbations c) a system of equations for describing the growth of these 
perturbations, their linear approximations and analytical solutions d) statistical 
quantities to describe the various uncertainties and selection effects e) nonlinear 
approximations of various kinds to obtain physical insight and f) full scale N 
body simulations in an expanding background. 

The following chapters detail the approaches we have taken and progress 
made, using the theoretical model sketched herein. 



Chapter 2 

Nonlinear gravitational 
clustering: Dreams of a paradigm 

Humor is the only test of gravity and gravity of humor - Aristotle 

2.1 Introduction 

As is clear from the first chapter the evolution of large number of particles under 
their mutual gravitational influence is a well-defined mathematical problem. But, 
if such a system occupies a finite region of phase space at an initial instant, and 
evolves via Newtonian gravity, then it does not reach any sensible 'equilibrium' 
state. The core region of the system will keep on shrinking and will be eventually 
be dominated by a few 'hard binaries'. Rest of the particles will evaporate away 
to large distances, gaining kinetic energy from the shrinking core [for a discussion 
of such systems, see [35] ] . 

The situation is drastically different in the presence of an expanding back- 
ground universe characterised by an expansion factor a(t). Firstly, the expansion 
tends to keep particles apart thereby exerting a civilising influence against New- 
tonian attraction. Secondly, it is now possible to consider an infinite region of 
space filled with particles. The average density of particles will contribute to the 
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expansion of the background universe and the deviations from the uniformity will 
lead to clustering. Particles evaporating from a local overdense cluster cannot 
escape to "large distances" but necessarily will encounter other deep potential 
wells. Naively, one would expect the local overdense regions to eventually form 
gravitationally bound objects, with a hotter distribution of particles hovering 
uniformly all over. As the background expands, the velocity dispersion of the 
second component will keep decreasing and they will be captured by the deeper 
potential wells. Meanwhile, the clustered component will also evolve dynamically 
and participate in, e.g mergers. If the background expansion and the initial con- 
ditions have no length scale, then it is likely that the clustering will continue in 
a hierarchical manner ad infinitum. 

Most of the practising cosmologists will broadly agree with the above picture 
of gravitational clustering in an expanding universe. It is, however, not easy to 
translate these concepts into a well-defined mathematical formalism and provide 
a more quantitative description of the gravitational clustering. One of the key 
questions regarding this system which needs to be addressed is the following: 
Can one make any general statements about the very late stage evolution of the 
clustering ? For example, does the power spectrum at late times 'remember' the 
initial power spectrum or does it possess some universal characteristics which are 
reasonably independent of initial conditions ? [This question is closely related to 
the issue of whether gravitational clustering leads to density profiles which are 
universal. [32]]. 

In this chapter we address some of the above issues and show that it is 
possible to provide (at least partial) answers to these questions based on a simple 
paradigm. The key assumption we shall make is the following: Let ratio between 
mean relative pair velocity v(a,x) and the negative hubble velocity (—ax) be 



denoted by h(a, x) and let £(a, x) be the mean correlation function averaged over 
a sphere of radius x. We shall assume that h(a, x) depends on a and x only through 
£(a, that is, h(a,x) = h[£(a,x)]. With such a minimal assumption, we will 
be able to obtain several conclusions regarding the evolution of power spectrum 
in the universe. Such an assumption was originally introduced — in a different 
form — by Hamilton [19]. The present form, as well as its theoretical implications 
were discussed in [34], and a theoretical model for the scaling was attempted by 
Padmanabhan [37]. It must be noted that simulations indicate a dependence of 
the relation h(a,x) = h[£(a,x)] on the intial spectrum and also on cosmological 
parameters ( [43] , [44] , [42] , [30]). Most of our discussion is independent of this fact 
or can be easily generalised to such cases. When we need to use an explicit form 
for h we shall use the original ones suggested by Hamilton [19] because of its 
simplicity. 

Since this chapter addresses several independent but related questions, we 
provide here a brief summary of how it is organised. Section 2.1 studies some as- 
pects of nonlinear evolution based on the assumption mentioned above. We begin 
by summarising some previously known results in subsection 2.2.1 to set up the 
notation and collect together in one place the formulas we need later. Subsection 
2.2.2 makes a brief comment about the critical indices in gravitational dynamics 
so as to motivate later discussion. In section 2.3, we discuss the relation between 
density profiles of halos and correlation functions and derive the conditions un- 
der which one may expect universal density profiles in gravitational clustering. 
In section 2.4 we show that gravitational clustering does not admit self similar 
evolution except in a very special case. We also discuss the conditions for ap- 
proximate self-similarity to hold. Section 2.5 discusses the question as whether 
one can expect to find power spectra which evolve preserving their shape, even in 
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the nonlinear regime. We first show, based on the results of section 2.4, that such 
exact solutions cannot exist. We then discuss the conditions for the existence of 
some approximate solutions. We obtain one prototype approximate solution and 
discuss its properties. The solution also allows us to understand the connection 
between statistical mechanics of gravitating systems in the small scale and evolu- 
tion of correlation functions on the large scale. Finally, section 2.6 discusses the 
results. 

2.2 General features of nonlinear evolution 

Consider the evolution of the system starting from a gaussian initial fluctua- 
tions with an initial power spectrum, P in (k). The fourier transform of the power 
spectrum defines the correlation function £ (a, x) where a oc t 2 / 3 is the expansion 
factor in a universe with Q = 1. It is more convenient to work with the average 
correlation function inside a sphere of radius x, defined by 

e(a,x) = -| [ X ti(a,y)y 2 dy (2.1) 
x a Jo 

This quantity is related to the power spectrum P(a, k) by 
_ 3 r°° dk 

£(x : a) = — 5—5-/ —P(a, k) [sm(k x) — k x cos(/c x)] (2.2) 
2;, .r 1 Jo k 

with the inverse relation 

P{a,k) = — / dx x £(a, x) [sinffc x) — k x cos(/c x)} (2.3) 
ok Jo 

In the linear regime we have £l(o, x) oc a 2 ^ in (ai,x). 

We now recall that the conservation of pairs of particles gives an exact equa- 
tion satisfied by the correlation function [45]: 
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where v(a,x) denotes the mean relative velocity of pairs at separation x and 
epoch a. Using the mean correlation function £ and a dimensionless pair velocity 
h(a,x) = —(v/ax), equation (2.4) can be written as 

<aL- A aL»< 1 + «»- 3A < 1 + «» < 2 - 5 > 

This equation can be simplified by introducing the variables 

A = In a, X = \nx, D(X, A) = ln(l + £) (2.6) 
in terms of which we have [34] 

3D 3D 

M -h(A,X)— = 3h(A,X) (2.7) 

At this stage we shall introduce our key assumption, viz. that h depends on 
(A, X) only through £ (or, equivalently, D). Given this single assumption, several 
results follow which we shall now summarise. 

2.2.1 Formal solution 

Given that h = /i[£(a, x)}, one can easily integrate the equation (2.5) to find 
the general solution [34]. The characteristics of this equation (2.5) satisfy the 
condition 

x 3 (l + £) = / 3 (2.8) 

where I is another length scale. When the evolution is linear at all the relevant 
scales, ( C 1 and I ~ x. As clustering develops, £ increases and x becomes 
considerably smaller than I. The behaviour of clustering at some scale x is then 
determined by the original linear power spectrum at the scale / through the "flow 
of information" along the characteristics. This suggests that we can express the 
true correlation function £(a, x) in terms of the linear correlation function £l(o, /) 
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evaluated at a different point. This is indeed true and the general solution can 
be expressed as a nonlinear scaling relation (NSR, for short) between £l(o, I) and 
£(a, #) with / and x related by equation(2.8). To express this solution we define 
two functions V(z) and U(z) where V(z) is related to the function h(z) by 

and U(z) is the inverse function of V(z). Then the solution to the equation (2.5) 
can be written in either of two equivalent forms as: 

fax) =u[t L (a,l)] ; £ L (a,l) = V [fax)] (2.10) 

where I 3 = x 3 (l + £) [34]. Given the form of h(£) this allows one to relate the 
nonlinear correlation function to the linear one. 

From general theoretical considerations [37] it can be shown that V(z) has 
the form: 

'1 (z«l) 
V(z) = \ z 1 / 3 (1 < z < 200) (2.11) 
, z 2 ' 3 (200 < z) 

In these three regions h(z) m [(2z/3), 2, 1] respectively. We shall call these 
regimes, linear, intermediate and nonlinear respectively. More exact fitting func- 
tions to V(z) and U(z) have been suggested in literature. ([19], [30], [43]). When 
needed in this chapter, we shall use the one given in Hamilton [19]: 

/ 1 + 0.0158 * 2 + 0.000115 * 3 \ 1/3 
{Z> Z \ 1 + 0.926 z 2 - 0.0743 z 3 + 0.0156 z A J l ' ' 

z + 0.358 z 3 + 0.0236 z 6 
U[Z) ~ 1 + 0.0134*3 + 0.0020 z 9 / 2 { ' 

Equations (2.10) and (2.12,2.13) implicitly determine £(a,x) in terms of £l(ci,x). 
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2.2.2 Critical indices 

These NSR already allow one to obtain some general conclusions regarding the 
evolution. To do this most effectively, let us define a local index for rate of 
clustering by 

..M)^ (2.14) 
o In a 

which measures how fast (,(a,x) is growing. When £(a,x) <C 1, then n a = 2 
irrespective of the spatial variation of £ (a, x) and the evolution preserves the 
shape of £(a, x). However, as clustering develops, the growth rate will depend on 
the spatial variation of £(a, x). Defining the effective spatial slope by 

- K// ( , l)+ 3] = ^f^> (2.15) 

one can rewrite the equation (2.5) as 

n a = M^-n e// ) (2.16) 

At any given scale of nonlinearity, decided by £(a, x), there exists a critical spatial 
slope n c such that n a > 2 for n e // < n c [implying rate of growth is faster than 
predicted by linear theory] and n a < 2 for n e // > n c [with the rate of growth 
being slower than predicted by linear theory] . The critical index is fixed by setting 
n a = 2 in equation (2.16) at any instant. This feature will tend to "straighten out" 
correlation functions towards the critical slope. [We are assuming that £(a, x) has 
a slope that is decreasing with scale, which is true for any physically interesting 
case]. From the fitting function it is easy to see that in the range 1 & £ ^ 200, 
the critical index is n c ~ —1 and for 200 & £, the critical index is n c ~ —2 [10]. 
This clearly suggests that the local effect of evolution is to drive the correlation 
function to have a shape with (1/x) behaviour at nonlinear regime and (l/x 2 ) in 
the intermediate regime. Such a correlation function will have n a ~ 2 and hence 
will grow at a rate close to a 2 . 



2.3 Correlation functions, density profiles and 
stable clustering 

Now that we have a NSR giving £ (a, x) in terms of Z) we can ask the question: 
How does £(a, x) behave at highly nonlinear scales or, equivalently, at any given 
scale at large a ? 

To begin with, it is easy to see that we must have v = —ax or h = 1 for 
sufficiently large £(a,x) if we assume that the evolution gets frozen in proper 
coordinates at highly nonlinear scales. Integrating equation (2.5) with h — 1, we 
get £(a, x) = a 3 F(ax); we shall call this phenomenon "stable clustering". There 
are two points which need to be emphasised about stable clustering: 

(1) At present, there exists some evidence from simulations 

[42] that stable clustering does not occur in a Q = 1 model. In a formal sense, 
numerical simulations cannot disprove [or even prove, strictly speaking] the occur- 
rence of stable clustering, because of the finite dynamic range of any simulation. 

(2) . Theoretically speaking, the "naturalness" of stable clustering is often 
overstated. The usual argument is based on the assumption that at very small 
scales — corresponding to high nonlinearities — the structures are "expected to 
be" frozen at the proper coordinates. However, this argument does not take into 
account the fact that mergers are not negligible at any scale in an Q = 1 universe. 
In fact, stable clustering is more likely to be valid in models with Q < 1 — a 
claim which seems to be again supported by simulations [42]. 

//stable clustering is valid, then the late time behaviour of £(a, #) cannot 
be independent of initial conditions. In other words the two requirements: (i) 
validity of stable clustering at highly nonlinear scales and (ii) the independence 
of late time behaviour from initial conditions, are mutually exclusive. This is 
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most easily seen for initial power spectra which are scale-free. If Pi n (k) oc k n so 
that £l((i,x) oc a 2 x~( n+3 \ then it is easy to show that £(a,x) at small scales will 
vary as 

_ 6 3(n+3) _ 

£(a,x) oc a—x~^+^; (£ > 200) (2.17) 

if stable clustering is true. Clearly, the power law index in the nonlinear regime 
"remembers" the initial index. The same result holds for more general initial 
conditions. 

What does this result imply for the profiles of individual halos? To answer 
this question, let us start with the simple assumption that the density field p(a, x) 
at late stages can be expressed as a superposition of several halos, each with some 
density profile; that is, we take 

p(a,x) = ^/(x-x,,a) (2.18) 

i 

where the i-th halo is centered at Xj and contributes an amount /(x — x i? a) at the 
location Xj [We can easily generalise this equation to the situation in which there 
are halos with different properties, like core radius, mass etc by summing over 
the number density of objects with particular properties; we shall not bother to 
do this. At the other extreme, the exact description merely corresponds to taking 
the /'s to be Dirac delta functions]. The power spectrum for the density contrast, 
S(a,x) = (p/pb — 1), corresponding to the p(a, x) in (2.18) can be expressed as 



P(k,a) oc (a 3 |/(k, 



a 



exp — ik ■ Xj(a) 



(2.19) 



oc (a 3 |/(k,a)|) P cent (k,a) (2.20) 

where P cent (k, a) denotes the power spectrum of the distribution of centers of the 
halos. 

If stable clustering is valid, then the density profiles of halos are frozen in 
proper coordinates and we will have /(x — x,, a) = /(a(x — x,)); hence the fourier 
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transform will have the form /(k, a) = f(k/a). On the other hand, the power 
spectrum at scales which participate in stable clustering must satisfy P(k, a) = 
P(k/a) [This is merely the requirement £(a, x) = a 3 F(ax) re-expressed in fourier 
space]. From equation (2.20) it follows that we must have P ccnt (k, a) = constant 
independent of k and a at small length scales. This can arise in the special 
case of random distribution of centers or — more importantly — because we are 
essentially probing the interior of a single halo at sufficiently small scales. [Note 
that we must necessarily have P ccn t ~ constant, for length scales smaller than 
typical halo size, by definition]. We can relate the halo profile to the correlation 
function using (2.20). In particular, if the halo profile is a power law with / oc r _e , 
it follows that the £(a,x) scales as x^ 1 ( [29], [53]) where 

7 = 2e-3 (2.21) 

Now if the correlation junction scales as [— 3(n + 3)/(n + 5)], then we see that 
the halo density profiles should be related to the initial power law index through 
the relation 

e = 3 -^> (2.22) 
n + 5 v ' 

So clearly, the halos of highly virialised systems still "remember" the initial power 

spectrum. 

Alternatively, one can try to "reason out" the profiles of the individual halos 
and use it to obtain the scaling relation for correlation functions. One of the 
favourite arguments used by cosmologists to obtain such a "reasonable" halo 
profile is based on spherical, scale invariant, collapse. It turns out that one 
can provide a series of arguments, based on spherical collapse, to show that - 
under certain circumstances — the density profiles at the nonlinear end scale as 
[— 3(n + 3)/(n + 5)]. The simplest variant of this argument runs as follows: If we 
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start with an initial density profile which is r~ a , then scale invariant spherical 
collapse will lead to a profile which goes as r~@ with (3 = 3a/ (1 + a) [37]. Taking 
the intial slope as a = (n + 3)/2 will immediately give (5 = 3(n + 3)/(n + 5). [Our 
definition of the stable clustering in the last section is based on the scaling of the 
correlation function and gave the slope of [— 3(n + 3)/(n + 5)] for the correlation 
function. The spherical collapse gives the same slope for halo profiles.} In this 
case, when the halos have the slope of e = 3(n + 3)/(n + 5), then the correlation 
function should have slope 

3(n+l) 

Once again, the final state "remembers" the initial index n. 

Is this conclusion true? Unfortunately, simulations do not have sufficient 
dynamic range to provide a clear answer but there are some claims [32] that the 
halo profiles are "universal" and independent of initial conditions. The theoretical 
arguments given above are also not very rigourous (in spite of the popularity they 
seem to enjoy!). The argument for correlation function to scale as [— 3(n+3)/(n+ 
5)] is based on the assumption of h — 1 asymptotically, which may not be true. 
The argument, leading to density profiles scaling as [— 3(n+3)/(n+5)], is based on 
scale invariant spherical collapse which does not do justice to nonradial motions. 
Just to illustrate the situations in which one may obtain final configurations which 
are independent of initial index n, we shall discuss two possibilities: 

(i) As a first example we will try to see when the slope of the correlation 
function is universal and obtain the slope of halos in the nonlinear limit using 
our relation (2.21). Such an interesting situation can develop if we assume that 
h reaches a constant value asymptotically which is not necessarily unity. In that 
case, we can integrate our equation (2.5) to get £(a, x) = a 3h F[a h x] where h now 
denotes the constant asymptotic value of of the function. For an initial spectrum 
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which is scale-free power law with index n, this result translates to 

£(a, x) oc a^x' 1 (2.24) 

where 7 is given by 

3h(n + 3) 

^srkri) (2 - 25) 

We now notice that one can obtain a 7 which is independent of initial power law 
index provided h satisfies the condition h(n + 3) = c, a constant. In this case, 
the nonlinear correlation function will be given by 

£(a,x) oc a( 2 + c H™+ a )ar^ (2.26) 

The halo index will be independent of n and will be given by 

(2,7) 

Note that we are now demanding the asymptotic value of h to explicitly depend on 
the initial conditions though the spatial dependence of £(a,x) does not. In other 
words, the velocity distribution — which is related to h — still "remembers" the 
initial conditions. This is indirectly reflected in the fact that the growth of £ (a, x) 
- represented by a 6c /^ 2+c ^ n+3 ^ — does depend on the index n. 

As an example of the power of such a — seemingly simple — analysis note the 
following: Since c > 0, it follows that e > (3/2); invariant profiles with shallower 
indices (for e.g with e = 1) are not consistent with the evolution described above. 

(ii) For our second example, we shall make an ansatz for the halo profile and 
use it to determine the correlation function. We assume, based on small scale 
dynamics, that the density profiles of individual halos should resemble that of 
isothermal spheres, with e = 2, irrespective of initial conditions. Converting this 
halo profile to correlation function in the nonlinear regime is straightforward and 
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is based on equation (2.21): If e = 2, we must have 7 = 2e — 3 = 1 at small 
scales; that is £(a, x) oc x~ l at the nonlinear regime. Note that this corresponds 
to the critical index at the nonlinear end, n e ff = n c = — 2 for which the growth 
rate is a 2 — same as in linear theory. [This is, however, possible for initial power 
law spectra, only if e = 1, i.e h(n + 3) = 1 at very nonlinear scales. Testing the 
conjecture that h(n + 3) is a constant is probably a little easier than looking for 
invariant profiles in the simulations but the results are still uncertain]. 

The corresponding analysis for the intermediate regime, with 1 ^ £(a,x) & 200, 
is more involved. This is clearly seen in equation (2.20) which shows that the 
power spectrum [and hence the correlation function] depends both on the fourier 
transform of the halo profiles as well as the power spectrum of the distribution 
of halo centres. In general, both quantities will evolve with time and we cannot 
ignore the effect of P ccnt (k,a) and relate P{k,a) to f(k,a). The density profile 
around a local maxima will scale approximately as p oc £ while the density profile 
around a randomly chosen point will scale as p oc £ x / 2 . [The relation 7 = 2e — 3 
expresses the latter scaling of £ oc p 2 ] . There is, however, reason to believe that 
the intermediate regime (with 1 ^ £ ^ 200) is dominated by the collapse of 
high peaks [37] . In that case, we expect the correlation function and the density 
profile to have the same slope in the intermediate regime with £(a,x) oc (1/x 2 ). 
Remarkably enough, this corresponds to the critical index n e ff = n c = — 1 for 
the intermediate regime for which the growth is proportional to a 2 . 

We thus see that if: (i) the individual halos are isothermal spheres with (l/x 2 ) 
profile and (ii) if £ oc p in the intermediate regime and £ oc p 2 in the nonlinear 
regime, we end up with a correlation function which grows as a 2 at all scales. 
Such an evolution, of course, preserves the shape and is a good candidate for the 
late stage evolution of the clustering. 
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While the above arguments are suggestive, they are far from conclusive. It 
is, however, clear from the above analysis and it is not easy to provide unique 
theoretical reasoning regarding the shapes of the halos. The situation gets more 
complicated if we include the fact that all halos will not all have the same mass, 
core radius etc and we have to modify our equations by integrating over the 
abundance of halos with a given value of mass, core radius etc. This brings in 
more ambiguities and depending on the assumptions we make for each of these 
components [e.g, abundance for halos of a particular mass could be based on 
Press-Schecter or Peaks formalism], and the final results have no real significance. 
It is, therefore, better [and probably easier] to attack the question based on 
the evolution equation for the correlation function rather than from "physical" 
arguments for density profiles as done next. 

2.4 Self-similar evolution 

Since the above discussion motivates us to look for correlation functions of the 
form £(a,x) = a 2 L(x), we will start by asking a more general question: Does 
equation (2.5) possess self-similar solutions of the form 

Z(a,x)=a?F(±;) = af l F(q) (2.28) 
a" 

where q = xa~ a ?. Defining Q = \nq = X — a A and changing independent 
variables to from (A, X) to (A,Q) we can tranform our equation (2.5) to the 
form: 

Using the relations (<9£/cM)q = (<9£/<9Q)a = (|/ F)(dF/dQ) we can rewrite 
this equation as 

«-3(l + fM) 1V K{Q) (2 .30) 



a + h(g) f FdQ 
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The right hand side of this equation depends only on Q and hence will vanish if 
differentiated with respect to A at constant Q. Imposing this condition on the 
left hand side and noticing that it is a function of £ (a, x) we get 

( || J i(Left Hand Side) = (2.31) 

To satisfy this condition we either need (i) (d£/dA)Q = /3£ = implying /3 = 
or (ii) the left hand side must be a constant. Let us consider the two cases 
separately. 

(i) The simpler case corresponds to (5 = which implies that £(a, x) = F(Q). 
Setting (3 = in equation (2.30) we get 



<$\_ 3 ( i + «Mf) (2 32) 



which can be integrated in a straightforward manner to give a relation between 
q = exp Q and £: 

= g (i + £)- 1/3 v(O^ 2 

Given the form of h[£(a, x)], this equation can be in principle inverted to deter- 
mine £ as a function of q = xa~ a . 

To understand when such a solution will exist, we should look at the limit of 
(<C 1. In this limit, when linear theory is valid, we know that h (2/3)£ [45]. 
Using this in equation (2.33) we get the solution to be ln£ = — (2/a) lng or 

£ oc q-i oc x~ia 2 oc a 2 x~ (n+3) (2.33) 

with the definition a = 2/(n + 3). This clearly shows that our solution is valid, 
if and only if the linear correlation function is a scale-free power law. In this 
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case, of course, it is well known that solutions of the type £(a,x) = F(q) with 

2 

q = xa (™+ 3 ) exists. [Equation (2.33) gives the explicit form of the function F{q)\. 
This result shows that this is the only possibility. It should be noted that, even 
though we have no explicit length scale in the problem, the function £(q) - 
determined by the above equation — does exhibit different behaviour at different 
scales of nonlinearity. Roughly speaking, the three regimes in equation (2.11) 
translates into nonlinear density contrasts in the ranges S < 1, 1 < S < 200 and 
5 > 200 and the function £(g) has different characteristics in these three regimes. 
This shows that gravity can intrinsically select out a density contrast of 5 ~ 200 
which, of course, is well-known from the study of spherical tophat collapse. 

(ii) Let us next consider the second possibility, viz. that the left hand side 
of equation (2.30) is a constant. If the constant is denoted by //, then we get 
F = F and 

PZ-3(l + 0h(0= f iaZ + f iht (2.34) 
which can be rearranged to give 

h = iP ' afM) l (2.35) 
3 + (/i + 3)£ 1 ' 

This relation shows that solutions of the form £(a,x) = a 13 F(x/a a ) with (3^0 
is possible only if h[£(a, x)] has a very specific form given by (2.35). In this 
form, h is a monotonically increasing function of £(a,x). There is, however, 
firm theoretical and numerical evidence [19], [37] to suggest that h increases with 
£(a,x) first, reaches a maximum and then decreases. In other words, the h for 
actual gravitational clustering is not in the form suggested by equation (2.35). 
We, therefore, conclude that solutions of the form in equation (2.28) with [3 ^ 
cannot exist in gravitational clustering. 

By a similar analysis, we can prove a stronger result: There are no solutions 



of the form £(a, x) = £(x/F(a)) except when F(a) oc a a . So self-similar evolution 
in clustering is a very special situation. 

This result, incidentally, has an important implication. It shows that power- 
law initial conditions are very special in gravitational clustering and may not 
represent generic behaviour. This is because, for power laws, we have a strong 
constraint that the correlations etc can only depend on q = xa^ 2 ^ n+3 \ For more 
realistic — non-power law — initial conditions the shape can be distorted in a 
generic way during evolution. 

All the discussion so far was related to finding exact scaling solutions. It 
is however possible to find approximate scaling solutions which are of practical 
interest. Note that we normally expect constants like a, j3, \i etc to be of order 
unity while £(a, #) can take arbitrarily large values. If £(a,x) ^> 1 then equation 
(2.35) shows that h is approximately a constant with h = {(5 — aji)/ (/i + 3). In 
this case 

f(a, x) = a^F{q) oc aV oc a^'^x" oc a fe ^ +3 V (2.36) 

which has the form £(a, x) = a 3h F(a h x) which was obtained earlier by directly 
integrating equation (2.5) with constant h. We shall say more about such ap- 
proximate solutions in the next section. 

2.5 Units of the nonlinear universe 

Having reached the conclusion that exact solutions of the form £(a,x) = a 2 G(x) 
are not possible, we will ask the question: Are there such approximate solutions? 
And if so, how do they look like? We will see that such profiles — which we 
shall call "pseudo-linear profiles" — that evolve very close to the the above form 
indeed exist. In order to obtain such a solution and check its validity, it is better 
to use the results of section 2.2.1 and proceed as follows: 



We are trying to find an approximate solution of the form £(a, x) = a 2 G(x) 
to equation (2.5). Since the linear correlation function £l(g, %) does grow as a 2 at 
fixed x, continuity demands that £(a,x) = ^i(a,x) for all a and x. [This can be 
proved more formally as follows: Let £ = a 2 G(x) and £l = a 2 Gi(:r) for some range 
x\ < x < X2- Consider a sufficiently early epoch a = at at which all the scales 
in the range (x 1 ,x 2 ) are described by linear theory so that £(aj,x) = £l(o>ux)- 
It follows that Gi(x) = G(x) for all x\ < x < x 2 . Hence £(a,x) = £ L (a,x) for 
all a in X\ < x < x 2 . By choosing sufficiently small, we can cover any range 
(x±, x 2 ). So £ = £l for any arbitrary range. QED}. Since we have a formal relation 
(2.10) between nonlinear and linear correlation functions, we should be able to 
determine the form of G(x). 

To do this we shall invert the form of the linear correlation function £z,(a, /) = 
a 2 G(l) and write / = G^ 1 (a~ 2 ^ L ) = F(a~ 2 ^ L ) where F is the inverse function of 
G. We also know that the linear correlation function £ L (a,l) at scale I can be 
expressed as V[£(a, x)] in terms of the true correlation function £(a,x) at scale x 
where 

I = x(l + £(a,x)) 1/3 (2.37) 

So we can write 

I = p \^lR 

a 2 

But x can be expressed as x = F[£ L (a, x)/a 2 ]; Substituting this in (2.37) we have 

l + e] 1/3 (2.39) 

From our assumption £t,(a, x) = £(a,x) ; therefore this relation can also be 
written as 

a 2 



V[£(x,a)] 



(2.38) 



I = F 



£ L (a,x) 



(i+i) 



^V3 



(2.40) 
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Equating the expressions for / in (2.38) and (2.40) we get an implicit functional 
equation for F: 

-.1/3 





= F 


T 


. ° 2 . 







(1 + 5)' 



(2.41) 



which can be rewritten as 



F 


V(0/a 2 ; 


F 







;i + o 



1/3 



(2.42) 



This equation should be satisfied by the function F if we need to maintain the 
relation £,(a,x) = ^(0,2;). 

To see what this implies, note that the left hand side should not vary with a 
at fixed £. This is possible only if F is a power law: 



F(0 = AC 
which in turn constrains the form of V(£) to be 

v(0 = e(i + 1/3m 



(2.43) 



(2.44) 



Knowing the particular form for V we can compute the corresponding /i(£) from 
the relation 

(2.45) 



d£ 3 (1+0 HO 
For the V(£) considered in equation (2.44) we get 



h 



2£ 



(2.46) 



3 + (3 + l/m)i 

which is the same result obtained by putting f3 = 2 , a = in equation (2.28). 
We thus recover our old result — as we should — that exact solutions of the 
form ^(a, x) = £ L (a,x) = a 2 G(x) are not possible because the correct V(^) and 
h{£) do not have the forms in equations (2.44) and (2.46) respectively. But, as 
in the last section, we can look for approximate solutions. 
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We note from equation (2.44) that for £ 3> 1, we have 

V(0 = £( 1+1 /3 m ). F(e ") £l/m (2.47) 

This can be rewritten as 

V(0 = r; F(0 oc e 1 / 3 ^" 1 ); G(0 oc pA"" 1 ) (2.48) 

In other words if V(£) can be approximated as we have an approximate 
solution of the form 

f (a, a) = a 2 = a 2 a; 3 ^ 1 ) (2.49) 

Since the V in equation (2.12) is well approximated by the power laws in (2.11) 
so that 

V(0 oc ^ 3 (1 < l< 200) (2.50) 
oc f /3 (200 < (2.51) 

we can take v — 1/3 in the intermediate regime and z/ = 2/3 in the nonlinear 
regime. It follows from (2.48) that the approximate solution should have the form 

F{1) oc 4= (1 £ £ £ 20 °) (2-52) 

oc i (200 < (2.53) 

This gives the approximate form of a pseudo-linear profile which will grow as a 2 
at all scales. 

There is another way of looking at this solution which is probably more phys- 
ical and throws light on the scalings of pseudo- linear profiles. We recall that, in 
the study of finite gravitating systems made of point particles and interacting via 
Newtonian gravity, isothermal spheres play an important role. They can be shown 
to be the local maxima of entropy [35] and hence dynamical evolution drives the 
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system towards an (l/x 2 ) profile. Since one expects similar considerations to hold 
at small scales, during the late stages of evolution of the universe, we may hope 
that isothermal spheres with (l/x 2 ) profile may still play a role in the late stages 
of evolution of clustering in an expanding background. However, while converting 
the profile to correlation, we have to take note of the issues discussed in section 
2. In the intermediate regime, dominated by scale invariant radial collapse [37], 
the density will scale as the correlation function and we will have £ oc (l/x 2 ). 
On the other hand, in the nonlinear end, we have the relation 7 = 2e — 3 [see 
equation (2.21) ] which gives £ oc (l/x) for e = 2. Thus, if isothermal spheres are 
the generic contributors, then we expect the correlation function to vary as (l/x) 
and nonlinear scales, steepening to (l/x 2 ) at intermediate scales. Further, since 
isothermal spheres are local maxima of entropy, a configuration like this should 
remain undistorted for a long duration. This argument suggests that a £ which 
goes as (l/x) at small scales and (l/x 2 ) at intermediate scales is likely to be a 
candidate for pseudo-linear profile. It was found that this is indeed the case. 

To go from the scalings in two limits given by equation (2.52) to an actual 
profile, we can use some fitting function. By making the fitting function suffi- 
ciently complicated, we can make the pseudo-linear profile more exact. We shall, 
however, choose the simplest interpolation between the two limits and try the 
ansatz: 

F « = T^T^b) (2 - 54) 

where A and B are constants. Using the original definition I = F[£ L /a 2 ] and the 
condition that £ = £l, we get 
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This relation implicitly fixes our pseudo-linear profile. Solving for £, we get 

with L = 4A/_B 2 . Since this profile is not a pure power law, this will satisfy the 
equation (2.42) only approximately. We choose B such that the relation 

is satisfied to greatest accuracy at a = 1. This approximate profile works rea- 




Figure 2.1: The approximate solution to the functional equation determining the 
pseudo-linear profile. 



sonably well. Figures 2.1 and 2.2 show this result. In figure 2.1 we have plotted 
the ratio F(V(£) /a 2 )/ _F(£/a 2 ) on the x-axis and the function (1 + f) 1 / 3 on the 
y-axis. If the function in (2.56) satifies equation (2.42) exactly, we should get a 



45-degree line in the figure which is shown by a dashed line. The fact that our 
curve is pretty close to this line shows that the ansatz in (2.56) satisfies equation 
(2.42) fairly well. The optimum value of B chosen for this figure is B = 38.6. 
When a is varied from 1 to 10 3 , the percentage of error between the 45-degree 
line and our curve is less than about 20 percent in the worst case. It is clear that 
our profile in (2.56) satisfies equation (2.57) quite well for a dynamic range of 10 6 
in a 2 . 

Figure 2.2 shows this result more directly. We evolve the pseudo linear profile 
form a 2 = 1 to a 2 m 1000 using the NSR, and plot [£(a,x)/a 2 ] against x. The 
dot-dashed, dashed and two solid curves (upper one for a 2 = 100 and lower one 
for a 2 = 900) are for a 2 = 1,9,100 and 900 respectively. The overlap of the 
curves show that the profile does grow approximately as a 2 . Also shown are lines 
of slope —1 (dotted) and —2 (solid); clearly £ oc x~ l for small x and £ oc x~ 2 in 
the intermediate regime. 

We emphasis that we have chosen in equation (2.56) the simplest kind of 
ansatz combining the two regimes and we have used only two parameters A and 
B. It is quite possible to come up with more elaborate fitting functions which 
will solve our functional equation far more accurately but we have not bothered 
to do so for two reasons: (i) Firstly, the fitting functions in equation (2.11) 
for V(z) itself is approximate and is probably accurate only at 10-20 percent 
level. There has also been repeated claims in literature that these functions have 
weaker dependence on n which we have ignored for simplicity, (ii) Secondly, one 
must remember that only those £ which correspond to positive definite P(k) are 
physically meaningful. This happens to be the case our choice [which can be 
verified by explicit numerical integration with a cutoff at large x] but this may 
not be true for arbitrarily complicated fitting functions. Incidentally, another 



log(x/L) 



Figure 2.2: The dot-dashed, dashed and two solid curves (upper one for a 2 = 100 
and lower one for a 2 = 900) are for a 2 = 1,9, 100 and 900. The dotted straight line 
is of slope -1 and the solid one is of slope-2 showing both the l/x and 1/x 2 regions 
of the profile. 

simple fitting function for the pseudo-linear profile is 



with A' = B 2 and V = L/A. 

If a more accurate fitting is required, one can obtain it more directly from 
equation (2.16). Setting n a = 2 in that equation predicts the instantaneous 



£(a,x) 



= a 



2 



(2.58) 



(x/L')[(x/L') + l] 



spatial slope of £(a,x) to be 



<91n£(a, x) 



2 



3(1 + 



1 



(2.59) 



dlnx 



h[£(a,x)] 
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which can be integrated to give 

L Jt[ L ] £(2 -3h)-3h K J 

at a = 1 with L being an arbitratry integration constant. Numerical integration 
of this equation will give a profile which is varies as (1/x) at small scales and goes 
over to (1/x 2 ) and then to (1/x 3 ), (1/x 4 ).... etc with an asymptotic logarithmic 
dependence. In the regime £(a,x) > 1, this will give results reasonably close to 
our fitting function. 

It should be noted that equation (2.42) reduces to an identity for any F, in 
the limit £ — > since, in this limit V(z) z. This shows that we are free to 
modify our pseudo-linear profile at large scales into any other form [essentially 
determined by the input linear power spectrum] without affecting any of our 
conclusions. 

Finally, we will discuss a different way of thinking about pseudolinear profiles. 
In studying the evolution of the density contrast 5(a,x), it is conventional to 
expand in in term of the plane wave modes as 

5(a,x) = ]T5(a,k)exp(ik-x) (2.61) 

k 

In that case, the exact equation governing the evolution of 5 (a, k) is given by [45] 

da 2 2a da 2a? 

where A denotes the terms responsible for the nonlinear coupling between dif- 
ferent modes. The expansion in equation (2.61) is, of course, motivated by the 
fact that in the linear regime we can ignore A and each of the modes evolve 
independently. For the same reason, this expansion is not of much value in the 
highly nonlinear regime. 
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This prompts one to ask the question: Is it possible to choose some other set 
of basis functions Q(a, x), instead of exp ik ■ x, and expand 5(a, x) in the form 



so that the nonlinear effects are minimised ? Here a stands for a set of parameters 
describing the basis functions. This question is extremely difficult to answer, 
partly because it is ill-posed. To make any progress, we have to first give meaning 
to the concept of "minimising the effects of nonlinearity" . One possible approach 
we would like to suggest is the following: We know that when 5 (a, x) <C l,then 
<5(a,x) oc a -F(x) for any arbitrary -F(x); that is all power spectra grow as a 2 in 
the linear regime. In the intermediate and nonlinear regimes, no such general 
statement can be made. But it is conceivable that there exists certain special 
power spectra for which P(k, a) grows (at least approximately) as a 2 even in 
the nonlinear regime. For such a spectrum, the left hand side of (2.62) vanishes 
(approximately); hence the right hand side should also vanish. Clearly, such 
power spectra are affected least by nonlinear effects. Instead of looking for such 
a special P(k,a) we can, equivalently look for a particular form of £(a, x) which 
evolves as closely to the linear theory as possible. Such correlation functions 
and corresponding power spectra [which are the pseudo-linear profiles] must be 
capable of capturing most of the essence of nonlinear dynamics. In this sense, 
we can think of our pseudo-linear profiles as the basic building blocks of the 
nonlinear universe. The fact that the correlation function is closely related to 
isothermal spheres, indicates a connection between local gravitational dynamics 
and large scale gravitational clustering. 




(2.63) 



a 
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2.6 Results and Summary 

It seems reasonable to hope that the late stage evolution of collisionless point 
particles, interacting via Newtonian gravity in an expanding background, should 
be understandable in terms of a simple paradigm. This chapter tries to realise 
this within some well defined framework. It should be viewed as a first step in a 
new direction. 

There are three key points which emerge from this analysis. The first is 
the fact that we have been able to find approximate correlation functions which 
evolve preserving their shapes. We achieved this by looking at the structure 
of an exact equation which obeys certain nonlinear scaling relations. As we 
emphasised before, the existence of such special class of solutions to the equations 
of gravitational dynamics is an important feature. 

Secondly, we should take note of the role played by the "isothermal" profile 
(l/rr 2 ) in our solution. Such a profile can lead to correlation functions which 
go as (1/x) at small scales and (1/x 2 ) in the intermediate scales. If this profile 
is indeed "special" then one expects it to lead to a pseudo-linear profile for the 
correlation function. Our analysis shows that there is indeed good evidence for 
this feature. If one accepts this evidence, then the next level of enquiry would be 
to ask why (1/x 2 ) profiles are "special" . In the statistical mechanics of gravitating 
systems, one can show that these profiles arise as end stages of violent relaxation 
which operates at dynamical time scales. Whether a similar reasoning holds in 
an expanding background, independent of the index for power spectrum, is open 
to question. We emphasise that our equations, along with NSR, naturally lead 
to a pseudo-linear profile, which can be interpreted and understood in terms of 
isothermal density profiles for halos; we did not have to assume anything a priori 



regarding the halo profiles. 

In a more pragmatic way, one can understand the pseudo-linear profile from 
the dependence of the rate of growth of the correlation function on the local slope. 
The NSR suggest that £ grows (approximately) as a 6 / (ne - f - f+4 ) in the intermediate 
regime and as a 6//< " -e//+5 - ) in the nonlinear regime. This scaling shows that n e // = 

— 1 grows as a 2 in the intermediate regime and n e ff = —2 grows as a 2 in the 
nonlinear regime. This is precisely the form our pseudo-linear profile has. Also, 
in the intermediate regime, the correlation grows faster than a 2 if n e ff < — 1 and 
slower than a 2 if n e // > —1. The net effect is, of course, to straighten out a 
curved correlation and drive it to n — — 1. Similar effect drives the correlations 
to n — —2 in the nonlinear regime, [see [10] for a more detailed discussion of this 
aspect in the intermediate regime]. Of course, one still needs to understand the 
dependence of growth rate on the n e // from more physical considerations to get 
the complete picture. We have not addressed what is the timescale over which 
clustering can lead to the psuedo-linear profile assuming that it does. 

The last aspect has to do with what one can achieve using the pseudo-linear 
profiles. In principle, one would like to build the nonlinear density field through 
a superposition of pseudo-linear profiles but this is a mathematically complex 
problem. As a first step one should understand why the nonlinear term in equa- 
tion (2.62) is subdominant for such a profile. This itself is complicated since we 
have only fixed the power spectrum — but not the phases of the density modes 

— while the nonlinear terms do depend on the phase. 

High resolution numerical simulations serve to test the framework presented 
in this chapter but resource limitations hamper such experiments. The following 
chapter discusses the general theory of two dimensional simulations which may be 
used to explore these issues in greater detail by bypassing the resource constraints 
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plaguing 3-D simulations. 



Chapter 3 



A formal analysis of two 
dimensional gravity 

Computers make it easy to do a lot of things but most of the things they make it 
easier to do don't need to be done - Andy Rooncy 

3.1 Introduction 

The analytical explorations of the previous chapter has clearly shown that he 
equations describing the growth of density perturbations in the highly nonlinear 
stage are analytically intractable and hence large scale numerical simulations are 
resorted to for exploration of this regime. 

These N-Body simulations require large amount of computing resources (CPU, 
memory and storage space) if one is to get the requisite amount of dynamical 
range, i.e. good resolution in force and mass, large range in values of density 
etc. Time and resource constraints usually limit our ability to probe structure 
formation issues more deeply using computers, once the required resources are 
at the limits of technological feasibility. The key parameter which decides the 
feasibility level of numerical simulations is the size of a simulation, which — in 
turn — is characterised by: (i) The number of particles in the simulation volume, 
which is generally specified as N D , where D is the dimensionality of the simu- 
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lation (usually 2 or 3) and (ii) The number of mesh points (M) along any axis 
which determines the minimum length scale at which the results can be treated as 
reliable indicators of physical phenomena. In order to create a simulation volume 
that is a fair sample of the universe one needs about 10 7 particles and in order 
to have a high enough force resolution one needs to increase the number of grid 
points adequately [9]. Let us suppose we have 160 3 particles in three dimensions 
and our grid is 160 units on a side. Then, for the same amount of computational 
resources one can simulate a two dimensional situation with 2048 2 particles on 
a 2048 2 grid (160 3 « 2048 2 ). So, if we can extract useful (i.e. generalisable 
to the three dimensional case) physical insights from results in two dimensions, 
then simulations of two dimensional gravity will be helpful. This hope has led 
to a large number of two dimensional simulations in the field of gravitational 
clustering ( [7], [56], [50], [6], [51]). 

There are three ways in which two dimensional gravity can be operationally 
defined and corresponding numerical simulations undertaken: (i) Consider a sys- 
tem of point particles in a three dimensional (expanding) background with the 
force of interaction being given by Newton's law of gravitation (i.e F oc 1/r 2 ). 
The initial positions and velocities of the particles are such that they all lie in 
the same plane and all the velocities are in the plane i.e there are no velocity 
components orthogonal to the plane. This system will evolve with the particles 
being confined to the plane with clustering occurring in the plane. Thus, we have 
a two dimensional clustering scenario, (ii) Another system we can consider con- 
sists of infinite, thin 'needles' located parallel to each other. The mass elements 
in the 'needles' still interact through the 1/r 2 force, but the interaction between 
'needles' (obtained by summing over the mass elements) is given by a 1/r force. 
In this case as well, the background space expands uniformly in three dimensions. 



The two dimensional clustering that we study is the clustering of these 'needles', 
examined by taking a slice orthogonal to the 'needles', (iii) The third possibility 
involves writing down the Einstein's equations in two dimensions, finding the 
homogeneous and isotropic cosmological solution, taking the Newtonian limit (in 
which the potentials due to density perturbations and background metric can be 
superposed), finding the corresponding perfect fluid equations and solving them. 
In this case, we will have a background spacetime expanding in two dimensions 
unlike the other two cases. (There is yet another, fourth possibility, which can 
be defined only in an ad hoc manner discussed towards the end of the chapter.) 

The first case is not of much interest for cosmological simulations since the 
system is anisotropic, confined to a single plane and the clustering takes place 
in a specific plane only because the initial conditions were specifically selected to 
give this result. Hence we will not discuss it and it is mentioned here only for 
completeness. The way simulations in two dimensions are carried out usually is 
by simulating the second case and then defining the 'particles' as the intersection 
of the 'needles' with any plane orthogonal to them. In this case - as in the first 
case - the background spacetime expands in three dimensions (for a flat dust 
dominated universe the scale factor a(t) goes as t 2 ^ 3 ). The clustering that we 
observe and quantify in two dimensions, is basically the clustering of these needles 
in three dimensions. But this is also an anisotropic situation since the background 
spacetime is expanding in three dimensions. As an alternative we may try to write 
down the equations derived from Einstein's equations in two dimensions and 
examine how a system of particles interacting in a two dimensional expanding 
background spacetime is to be simulated. 

In the rest of this chapter we shall examine the third alternative. We will 
take a very general approach by developing the formal theory of (D + 1) gravity 
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and considering D = 2 as a special case [15]. 

The basic layout of the chapter is as follows: In section (3.2) we first define the 
analogue of Einstein's gravity in (D + 1) dimensions, discuss the Newtonian limit 
and its corresponding Poisson's equation and then go on to analyse the Friedmann 
metric in (D+l) dimensions for a flat universe with dust. In section (3.3) we write 
down the D dimensional fluid equations and obtain the equation governing the 
density perturbations. This equation is then solved in the linear approximation 
and using the Spherical Top Hat model. Then in section (3.4) and section (3.5) 
we specialise to the cases D = 3 and D = 2 respectively. Finally, in section (3.6) 
we summarise and discuss the implications of the results obtained in the earlier 
sections. 

3.2 Formal (D + 1) dimensional gravity 

We start our analysis of (D + l) dimensional (1 time dimension and D space 
dimensions) gravity from the action principle which we assume has the same 
form as that used in (3 + 1) dimensions. Using this action we construct the 
corresponding (D + 1) dimensional Einstein equations which will be subsequently 
used to study structure formation and spherical collapse. Thus, we begin with 
the action principle, 



where S g is the action for the gravitational field, <S m is the action for the matter 
fields, g is the determinant of the metric tensor g^, R is the Ricci scalar, k(D) 
is a suitable constant which can be, in general, a function of D (when D = 3, 
k = 8nG, G being the usual gravitational constant) and £ m is the lagrangian den- 
sity for the matter fields. The metric signature we adopt is (+, — ). 




(3.1) 
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We adopt the following convention regarding indices. Latin alphabets i,j,k... 
are used to represent (D + 1) dimensional indices which take on the values 
(0, 1, 2, . . . , D) while greek letters are used to denote D dimensional indices tak- 
ing on the values (1, 2, . . . , D). Varying the total action S with respect to g ik we 
obtain Einstein's equations, 

Gik = Rik — 2 9ik R = ~~j^T ik (3.2) 
where T ik is the energy momentum tensor of the matter fields and is defined by 

i [\~\ r p. = d(vfe-) _ a ( d(vfe \ 

2^ m lk dg ik dx l \d{dg ik /dx l )j 1 } 

Gik is the Einstein tensor and R ik is the usual Ricci tensor. Note that the (1/2) 
that appears in Einstein's equations arises due to the square root in the term 
\g\ and has nothing to do with the dimension of the spacetime. 
We will use the above equations in the subsections to follow. In subsec- 
tion (3.2.1), we will study the Newtonian limit of the metric tensor and then 
construct the corresponding Poisson equation that relates the Newtonian gravi- 
tational field (j) to the matter density p. Then, in subsection (3.2.2), we analyse 
the Friedmann metric in (D + 1) dimensions and the corresponding Newtonian 
limit of this metric is derived. 

3.2.1 Poisson equation in D dimensions 

In this section, we derive the Poisson equation relating the gravitational potential 
4> to the matter density p. We keep all factors of c since the Newtonian limit 
involves the limit c — > oo. The analysis here follows closely the treatment in [28]. 
Consider the metric. 

ds 2 = (l + ^) c 2 dt 2 - dl 2 (3.4) 
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where is a function of space and time with dimensions of velocity square. The 
term dl 2 is the D dimensional spatial line element given by the formula 

D 

dl 2 = Y,{dx a f (3.5) 

a=l 

We will now show that the metric written above is the Newtonian limit of Ein- 
stein's gravitational equations. We do this by showing that, in the Newtonian 
limit, the equation of motion of a particle follows Newton's force law with the 
force — mV0... In relativistic mechanics, the motion of a particle of mass m is 
determined by the action function S 



S = —mc j 'ds = —mc J cdt^ ^1 H — j -^j (3.6) 

where v 2 is the square of the magnitude of the particle's velocity in D dimen- 
sions. In arriving at the second equality we have used the form of the metric in 
equation (3.4). In the limit c — > oo, the action S can be approximated as 

S « -mc 2 Jdt(l + -^J—J = J dt {^-mc 2 + ^mv 2 - mcj^j (3.7) 

The equation of motion for the particle can be immediately written down and we 
obtain 

= -mV0 (3.8) 

where v is the velocity vector in D space dimensions. Thus, Newton's force 
law is recovered in the non-relativistic limit and from this we conclude that the 
metric given in equation (3.4) is the Newtonian limit of Einstein's gravitational 
equations with acting as the Newtonian gravitational potential. 

The relation between and the mass density p is found by taking the c — > oo 
limit of Einstein's equations. This procedure, in (3 + 1) gravity, determines the 
constant k(D) since the Poisson equation is explicitly known. In other dimensions 
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however, a definite criterion, like Gauss's law for example, must be imposed in or- 
der to determine k(D). We now consider the limit c — > oo of Einstein's equations 
in the following manner. First, we use the line element given in equation (3.4) to 
calculate the Ricci tensor component i?oo 

^oo = ^pnp^(M W) - ^ W ( 3 - 9 ) 

where the summation convention has been invoked in the above equation and the 
sum over p is only over the spatial dimensions. Then, using equation (3.2), we 
obtain, 

where we have used the fact that gikg lk — D + 1 and assumed D ^ 1. Thus, 
Einstein's equations can be written in the equivalent form, 

Rik = ^ 4 - {Tik — jj _ - gikTj . (3-11) 

where T is the trace of T ik . The energy momentum tensor of point particles is 
Tik = pc 2 UiUk where p is the mass density and Ui is the four velocity. Since, in the 
non-relativistic limit, the macroscopic motion is slow, the space components of 
Ui can be neglected and only the time component should be retained. Therefore, 

,2 



u o — t/9oo an d ~ for all \i. Consequently, only Too — doopc is non-zero. 



Substituting for T ik into equation (3.11) and using the expression in equation (3.9) 
for Rqq, we get, 



^-(§z|)( 1 + |)^-(§z|)^ 



That is, 



V 2 = (^-|) k(D) P (3.13) 

where V 2 is the usual Laplacian operator in D dimensions. This equation is the 
Poisson equation in D dimensions. Note that when substituting for the value 
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of Roo from equation (3.9), we neglected the first term in comparison with the 
second since the former is of order c -4 while the latter is only of order c~ 2 . 

3.2.2 Friedmann Universe in (D + 1) dimensions 

Let us next consider the maximally symmetric Robertson- Walker metric in (D+l) 
dimensions, specialising to flat space with k = (we set c = 1 in this and in 
subsequent sections), 

ds 2 = dt 2 -a 2 (t)dl 2 (3.14) 

where a(t) is the scale factor and dl 2 is the D dimensional line element given in 
equation (3.5). Calculating the components of the Einstein tensor, we obtain, 

D(D - l)d 2 



G 



oo 



2a 2 

G n = G 22 = ... = G DD = {l-D)ad + 

(l-j)(D-l)d 2 (3.15) 

where a stands for da(t)/dt and similarly d is the second derivative of a{t) with 
respect to time. All the other components are zero. For consistency, the energy 
momentum tensor must have the form T\ = diag(p, —p, —p, —p, . . .) where p is 
the matter density and p is the pressure. 

Substituting in Einstein's equations, we obtain, 

D(D -l)d 2 , . . . 

= k(D) P (3.16) 



2a 2 

a D-2d 2 k(D)p 
a + 2 ^ ~D-1 



+ ^- 2 = (3-17) 



The above two equations, together with the equation of state in the form p = p(p) 
completely specify the system. Solving these three equations, we can determine 
a(t), p(t) and subsequently pit). Combining equations (3.16,3.17), we get the 
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single equation, 

a k(D) 



We now specialise to the case of pressureless dust with the equation of state 
p = 0. Using the principle of conservation of energy and momentum expressed 
by the relation 

T% = (3.19) 

we derive the following relation, 



,9:r fc V vm 7 2dx i 

= 47^(« Dt ")-^ tM = ° ( 3 - 20 ) 
a D Q x k V V 2 dx l v ' 

Noting that the only non-zero component of T\ is T° = p we finally get 

paP = constant = C\ (3-21) 

Substituting the above relation into equation (3.16), and solving for a(t) and 
subsequently for p(t), we obtain the solutions, 

Let us next consider the Newtonian limit of the Friedmann metric. This limit is 
important because the length scales of interest in structure formation are small 
compared to the Hubble radius and the velocities in the system are also much 
smaller than c. This permits us to study the formation of large scale structures 
in the universe in a Newtonian framework where the effective potential due to 
the expanding background universe, $frW) and the potential due to the density 
perturbations, ip, can be simply superposed. In order to obtain $frw> we first 
recast the Friedmann metric in equation (3.14) into the more convenient form 

ds 2 = dt 2 - a 2 (t) (dX 2 + X 2 dVl 2 ) (3.23) 
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where X is the radial distance in D dimensions and is the corresponding solid 
angle. We then apply the transformations (see [38], pg 80,346) 



where only terms up to quadratic in X are retained. Direct calculations, correct 
upto this order, transforms the Friedmann line element to the form, 



which upon comparison with the metric in equation (3.4) gives the equivalent 
Newtonian potential $frw in D dimensions as 



We will now use the results developed in the last two subsections to study struc- 
ture formation and spherical collapse using the STH model. 

3.3 Structure formation in D dimensions 

Having determined the form of the Poisson equation in the Newtonian limit and 
analysed the Friedmann equations in (D + 1) dimensions, we proceed to derive 
the equation for the growth of inhomogeneities in the expanding universe. After 
this we consider a specific model, the STH model, to study spherical collapse of 
matter. 

3.3.1 Equation for density perturbations in D dimensions 

Let us assume that matter in the universe is a perfect, pressureless fluid with 
density p m and flow velocity U. We can formally write down the D dimensional 
fluid equations describing a perfect fluid in an external potential field <3> tot in a 




(3.24) 




(3.25) 




(3.26) 
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proper coordinate system labelled by the D dimensional vector r. Therefore, we 
have, 

/ 9 p n 



dt j 
fdXJ 



+ V r -( Pm U) = (3.27) 



l + (U-V r )U = -V r $ tot (3.28) 

V / r 

where equation (3.27) is the usual continuity equation while equation (3.28) is 
the Euler equation for the fluid. The potential in equation (3.28), <& tot , is the 
total external Newtonian potential 

$tot = $frw + V? (3-29) 

where $frw is the background potential associated with the smooth background 
matter density pbm and is given in equation (3.26) while ip is the potential caused 
by density perturbations (p m — pbm)- The potential if satisfies the Poisson equa- 
tion given in equation (3.13). Thus, 

V rV = (§rf) < D )(P™ ~ Pbm) = ) <D)p hm 5 (3.30) 



where 5 is the density contrast defined by 

r Pm — Pbm 



(3.31) 



Pbm 

We now transform to comoving coordinates defined by x = r/a(t) and define the 
peculiar velocity v by the relation 

U = H(t)r + v = ax + v (3.32) 

where v = ax and H(t) = (a /a). Then, equation (3.27) and equation (3.28) 
become 

(^^j + DHp m + ^V x • (p m v) = (3.33) 

fdv\ 1 1 

— + Hv+- v-V x )v = — V x ^ (3.34) 
V at I a a 
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where we have used equation (3.26) to substitute for $frw- Similarly, in co- 
moving co-ordinates, equation (3.30) reduces to 

= ^) a W (3.35) 

Using p m = Pbm(l + 8), transforming the time variable from t to a(t) and defining 
a new velocity variable u by 

u = — = — (3.36) 
da aa 

we can obtain equations for 6(a) and u(a). Therefore, using equation (3.21) 
and performing the transformations, equation (3.33) and equation (3.34) further 
reduce to, 

Of 

— + V x -[u(l + 5)] = (3.37) 

9 <9u ( d 2 \ 9/ , 1 . 

a 2 — + a + 2— u + a 2 u • V x u = — -V x ^ 3.38 
oa \ a J a 1 

Now, we use the Friedmann equations in equations (3.16,3.17) with p replaced by 
Pbm and with p = to substitute for a in the above equation. Further, we define 
a new potential ^ by the relation 

so that, upon using equation (3.35), one obtains, 

where all reference to k(-D) has disappeared. Hence the final system of equations 
we need to tackle are, 

^ + V x -[u(l + 5)] = (3.41) 

^ + (u • V x )u = - Q —^-A [V x tf + u] (3.42) 
oa la 
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where A is given by the relation 



2k(D) \ a 2 p bm (*) . „ ._ ( D(D-1) \ a 2 ^ 

■ Pbm = 777- S Pc = /m -» (3-43) 



\D(D — 1) J a 2 rDm p c (t) ' rc \ 2k(D) J a 2 
For the k = universe, we will set A = 1. 

To proceed further and determine the equation satisfied by 5, we decompose 
the term d a up (where up is the /3th covariant component of the vector u and d a 
is short for d/dx a ) as 

cU/3 = ^ + O*/? + -^M a, /3 = (1, 2, . . . , D) (3.44) 

where a a/ 3 is the traceless, symmetric shear tensor, VL a p is the antisymmetric 
rotation tensor, 9 is the (trace) expansion and 8 a p is the Kronecker delta symbol. 
Then, equation (3.41) and equation (3.42) are combined by taking the divergence 
of equation (3.42) and using the above decomposition of d a up to obtain a single 
equation for 5. Straightforward algebra gives, 



da 2 V 2a J da \ 2a 2 J v 7 \ D J(l + 5)\da y 

(3.45) 

where a 2 = a a/3 a al3 and fi 2 = (l/2)fl a/3 fl a ^ . This equation is the full non- 
linear equation for 5. Apart from the obvious nonlinear terms containing 5 2 and 
(dS/da) 2 , the term (1 + 5)(a 2 — 2f2 2 ), which is the contribution from the shear 
and rotation, is also non-linear. The non-linear terms in 5 in the above equation 
render the equation unsolvable in general. Ignoring these non-linear terms to a 
first approximation, we can get a linear equation for 5(a), 

Assuming a power law solution for delta in the form 5 oc a p , we get, 

p = ^ ± -V9D 2 - 2AD + 16 (3.47) 
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as the required values for p. Notice that 8 has a growing mode as well as a 
decaying mode in general. The above solutions hold for all values of D > 1 in 
the linear regime. 

Though the full non-linear equation is not solvable, by neglecting the contri- 
bution from the shear and rotation terms and by using a suitable ansatz for 5, 
the resulting non-linear equation can be solved. We proceed to do this within the 
framework of the STH model in the next section. 

3.3.2 The Spherical Top Hat (STH) model 

In the STH (spherical collapse) model, we assume spherical symmetry by neglect- 
ing the shear and rotation terms in the equation for 5. With this assumption the 
S equation can be exactly solved. 

Transforming equation (3.45) by changing the independent variable back to t, 
dropping the rotation and shear terms and using the Friedmann equations given 
in equations (3.16,3.17), we get, 

d 2 5 hdS (D + l\ 1 fd5\ 2 (D-2\ r/i , AQ . 

We now define a function R(t) by the relation 

1 + 6 = ^ = C D R°(t) Phm (3 ' 49) 
where Cd = 2tt d ^ 2 /(DT[D/2j) is the volume of a unit sphere in D dimensions 
introduced for later convenience and M is a constant. The expression for 5 above 
can be rewritten using the relation p^a = p Q a^ from equation (3.21): 



1 + 5 



M \aA D . a D 



V (3-50) 



CdPqo-E 

where po an d ao are the matter density and scale factor at some (arbitrarily 
chosen) "present" epoch to- Substituting equation (3.50) in equation (3.48), we 
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get an equation for the growth of R(t) as, 

d 2 R D-2 k(D) M 



(3.51) 



dt 2 D(D-l) C D R D - 1 

[As an aside we may note that if the universe contains matter or fields with 
equations of state other than p = 0, the equation for R(t) becomes 

d 2 R D-2 K (D) M k(D) „„ , „ , „ , 

lF = -^T)^^-W^T) iXD - 2)p + Dp) ^ R (3 - 52) 

where the term ((D — 2)p + Dp) rest comes from the smoothly distributed compo- 
nent with p 7^ 0.] 

From the form of the equation of motion of R(t) we can give the following 
interpretation. Since the entire system considered above is spherically symmet- 
ric, we interpret R as the radius of a D dimensional spherical region containing 
a mass M. The equation of motion of R determines the motion of the surface of 
this region. In general, a spherical overdense region will be expected to initially 
expand because of the expansion of the background universe till the excess grav- 
itational force due to the overdensity of enclosed matter stops the expansion and 
causes the region to collapse back on itself. We will discuss the cases D = 3 and 
D = 2 in the subsequent sections and determine the differences in the behaviour 
of the growth of inhomogeneities. 

3.4 Summary of standard results in three di- 
mensions 

When D = 3, all the standard equations are recovered. First the Poisson equation 
satisfied by the Newtonian gravitational potential given by equation(3.13) reduces 
to the standard form, 

V 2 = AnGp (3.53) 
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where we have defined G by relating it to re(3) by k(3) = 8nG. Similarly equa- 
tions (3.41) and equation (3.42) reduce to (with A = 1), 

^ + V x -[u(l + 5)] = (3.54) 
oa 

^ + (u • V x )u = - A [V x * + u] (3.55) 
aa la 



while the equation for \l> becomes, 



V 2 ^ = -. (3.56) 
a 

In a similar manner, the 5 equation reduces to, 

d 2 5 3 d5 3 r/1 4 / cW\ 2 . 9 ^_ 9 . . 

-r^H : 7 <J(l + <5) = — — +(l + <J)(o- 2 -2Q 2 ) (3.57) 

da 2 2arfa 2a 2 V ; 3(1 + 5) \da) ' V ' 

and the solutions to the linear perturbation equation which is obtained by drop- 
ping the nonlinear terms and the [a 2 — 2Q 2 ) term, are 

3 

5 oc a p , p=l, -- (3.58) 

which are well known. The STH model for D = 3 also reduces to the standard 
form 

d 2 R GM AnG , j . 

which is again well known [38] . Therefore, it is seen that the full D dimensional 
equations reduce to the correct equations in three dimensions. Now we will go 
on to discuss the important case of D = 2. 

3.5 Two dimensional gravity 

If we naively consider the limit D — > 2 in the D dimensional equations, assuming 
that k(D) is finite in this limit, we obtain the following results. First, the Poisson 
equation (3.13) reduces to 

V 2 = (3.60) 
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The above result shows that in two dimensions, the gravitational potential does 
not couple to the matter density p. In structure formation, this means that 
inhomogeneities cannot grow since the perturbed potential ip is not related to 5 
at all. The second interesting result is that the background Newtonian potential 
$frw vanishes. This occurs because, referring back to equation (3.18), a = 
for pressureless dust and hence the background potential is zero. Further, the S 
equation reduces to 

d 2 S 2 d5 3 / d5\ 2 ,„ , , „,s 

^ + aTa-WTS){^) +d + ^ 2 -^ 2 )- (3.61) 

Linearising the equation as before by dropping the (a 2 — 2Vt 2 ) and (dS/da) 2 terms 
we obtain 

— -— f3 62) 

da 2 a da 

The solutions to the linearised equation are 

5oca p , p = 0,-l. (3.63) 

Thus only a constant or the decaying mode is present. This is consistent with 
the result that the perturbed gravitational potential does not couple to 5. If one 
considers the STH model, it is easy to see that the growth equation for R(t) 
reduces to 

d 2 R 

— = -K(2) Prest R (3.64) 

For, p rest = 0, the solution to the above equation is just R(t) = Bit + Bi where 
Bi, B 2 are constants. This is to be expected since there is no gravitational force 
which can lead to clustering and as a consequence the radius simply grows with 
time just like the background universe. Thus, if k(D) is finite in the limit D — > 2, 
it is not possible to have gravitational clustering that can grow with time. 
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We can, however, try some alternative approaches to examine whether it 
is possible to have a consistent physical picture of growing structures for two 
dimensional gravity. 

One possibility is that instead of assuming k(D) to be finite, let us assume 
that the expression (k,(D)(D — 2)) remains finite when D — > 2. This finite value 
can be fixed, for example, by invoking Gauss's theorem in D dimensions. This 
gives 

. * (D — 1\ 2tx d I 2 G 

K(D) = {—2) fm (365) 

Thus, k(D) — > 00 when D — > 2, but the Poisson equation acquires the form 

V 2 = 2nGp (3.66) 

This is, of course, the same form which is obtained by applying Gauss's law 
in two dimensions. Hence, as in three dimensions, the gravitational potential is 
determined by the matter density and thus inhomogeneities can in principle grow. 
There are, however, difficulties with this approach. To begin with, the constant 
factor c 4 /(2k(D)) in the action S in equation(3.1) vanishes for D = 2. But 
this is not too serious a problem. The gravitational part of the action certainly 
vanishes but because only the variations about the action are of significance this 
difficulty can be ignored. But a more serious problem arises when the solutions 
to the Friedmann equation are considered. The solutions for a(t) and p(t) in D 
dimensions are given in equation (3.22). Using equation (3.65), these reduce to, 

a(t) ( d- d/2 g Ci \ 1/D /D . {D-mm t - 2 

a{t) -[(D-2)T[D/2]) 1 ' P^~ C ' a ~ Dn D/2 G 1 

(3.67) 

When D — > 2 then, a — > 00 and p — > irrespective of the dependence on t. This 
implies that one cannot solve the equations describing the growth of structure in 
a consistent and non-singular way. Hence, we conclude that it is not possible to 



have a theoretical formulation of two dimensional gravity as the Newtonian limit 
to Einstein's equations in two dimensions. 

An alternative that remains is to use the Newtonian fluid equations in D 
dimensions directly and rewrite them for an expanding background with an ar- 
bitrary scale factor a(t). Note that a(t) is not obtained from the Friedmann 
equations and is completely arbitrary. We can superpose the potentials for the 
background universe and the perturbations in this case as before. The further 
assumptions we need to make are (i) the potential of the background universe 
$b g is of the form 

$ bg = -l% 2 (3.68) 
and (ii) the Poisson equation is given by 

V 2 = K(D) Phm 5 (3.69) 

where «(£>) = 2tt d / 2 G/T[D/2]. This form of «(£>) is obtained from the use 
of Gauss's law in D dimensions. We also need to specify how the background 
density pbm depends on time. In analogy with the usual Friedmann equations, 
we will assume Pb m a D = C\ where C\ is a constant. This gives an equation for S 
with an arbitrary scale factor a(t) as 



2 



d 2 5 (aa + 2a 2 \d5 1 r/i / D + 1\ 1 ( d5Y ... 2 nr>2 . 



da 2 \ ad 2 

(3.70) 

The above equation can be solved in any dimension D if the form of a(t) is given. 
This gives us a non-singular way to analyse growth of structures in D dimensions, 
including the case D = 2. But in three dimensions we observe that the above 
equation does not correctly reduce to equation (3.57). We may obtain the correct 
equation in three dimensions by making an additional ansatz, namely, that 

a) =— (3 - 71) 
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With this, equation (3.70) reduces to 

d 2 5 {Q-D\d5 D (D + l\ 1 fd5\ 2 M ... 2 oo2 , 

(3.72) 

Notice that the above equation differs from the earlier equation for 5 in D di- 
mensions, equation (3.45), in that there is a factor of (D — 2) missing from the 
coefficient of the 5(1 + 5) term. Therefore, the above equation does correctly 
reduce to equation (3.57) since (D — 2) equals unity when D = 3. When D = 2, 
equation (3.72) gives 

d 2 5 2 d5 1 r/1 rx 3 / d5\ 2 . w 2 _ oX 

S? + aS;-^ (1 + 5) = 2(TT5)Uj +d+*)(- 2 -2f! 2 ) (3-73) 

On linearising this equation by dropping the (a 2 — 2Q 2 ) term and the nonlinear 
terms and solving it we get both a growing mode as well as a decaying mode for 
5. The solutions are 

5 oc a\ q= (-1 ± VE)/2 (3.74) 

(It is interesting to note that one of the power law exponents is the golden ratio). 
The Spherical Top Hat (STH) equation in this case turns out to be 

f, GM lR /o _ x 
R = — + 1¥ (3.75) 

where M is the constant mass inside a 'spherical' shell of radius R. This equation, 
unfortunately, has no simple analytic solution. 

While this procedure leads to nontrivial results, it has many ad hoc assump- 
tions and cannot be obtained by taking appropriate limits of Einstein's theory in 
a systematic manner. Consequently it cannot be applied to numerical investiga- 
tions of 2 dimensional gravity with the confidence that the results will have some 
implications for the three dimensional case. 
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3.6 Results and Summary 

In this chapter we have analysed the case of two dimensional gravitational clus- 
tering starting from a formulation of the D dimensional Einstein's equations and 
taking the proper limits. The system of equations thus arrived at for a (D + 1) 
dimensional universe has been shown to reduce to the correct equations in three 
dimensions. But when the D — > 2 limit of these equations is taken, we are forced 
to conclude that irrespective of the value of k(D), a consistent two dimensional 
gravity theory in a cosmological context that supports growth of structures can- 
not be constructed. 

If k(2) is assumed to be finite, we observe that the coefficient in Poisson's 
equation goes to zero thus decoupling the potential from the density. This implies 
that perturbations do not grow but decay in time due to the expansion of the 
background spacetime. The alternative which is obtained by using the expression 
for k(D) given by equation (3.65) gives rise to solutions for the scale factor which 
are singular and therefore unacceptable. 

We have discussed all the ways in which two dimensional gravity may be 
simulated including an ad hoc procedure without a strong foundation which can 
give non-singular results as far as structure formation scenarios in two dimensions 
are concerned. The results presented here leads us to conclude that the only way 
to do a numerical simulation of two dimensional gravity is to simulate infinite 
'needles' in a background spacetime expanding in three dimensions and consider 
the 'particles' in the system to be intersections of the 'needles' with any plane 
orthogonal to them which is the approach taken in the next chapter. 



Chapter 4 

Scaling Relations in Two 
Dimensions 

When in doubt use brute force - Ken Thompson 

4.1 Gravitational Clustering: Two vs Three Di- 
mensions 

A number of attempts have been made in recent years to understand the evolu- 
tion of constructs like the two point correlation function using certain non-linear 
scaling relations (NSR) as discussed in chapter 2 ([19]), [34], [37]). These studies 
have shown that the relation between the non-linear and the linearly extrapolated 
correlation functions is reasonably model independent. This relation divides the 
evolution of correlation function into three parts [8]: the linear regime, the in- 
termediate regime and the non-linear regime. The evolution in the intermediate 
regime can be understood in terms of radial collapse around density peaks [37], 
if it is assumed that the evolution of profiles of density peaks follows the same 
pattern as an isolated peak. It is customary to invoke the hypothesis of stable 
clustering [45] to model the non-linear regime. A large number of studies have 
examined clustering in this regime and the general consensus is that the stable 
clustering limit does not exist [42]. 



83 



However, the limited dynamic range of currently available 3-dimensional N- 
Body simulations poses serious difficulties in investigating this problem in greater 
detail. It was pointed out [39] that we can circumvent this problem by simulating 
a two dimensional system, wherein a much higher dynamic range can be achieved. 
For example, since 1 60 3 ~ 2048 2 , the computational requirements are the same 
for a 2D simulation with box size of 2048 and 3D simulation with box size of 160. 
Assuming that one can reliably use, say, half of box size as good dynamic range 
we have a dynamic range of factor 1000 in 2D against a factor of about 80 in 3D. 
This allows us to probe higher nonlinearities in 2D compared to 3D. As long as we 
stick to generic features (like the non-linear scaling relations investigated here) 
which are independent of dimension, 2D has a definite advantage over 3D. Higher 
dynamic range is the basic motivation for studying gravitational clustering in two 
dimensions. 

When we go from three to two dimensions, we have, as discussed in chapter 
3, two different ways of modeling the system: 

(i) We can consider two dimensional perturbations in a three dimensional 
expanding universe. Here we keep the force between particles to be 1/r 2 and 
assume that all the particles, and their velocities, are confined to a single plane 
at the initial instant. 

(ii) We can study perturbations that do not depend on one of the three coor- 
dinates, i.e., we start with a set of infinitely long straight "needles" all pointing 
along one axis. The force of interaction falls as 1/r. The evolution keeps the 
"needles" pointed in the same direction and we study the clustering in an or- 
thogonal plane. Particles in the N-Body simulation represent the intersection of 
these "needles" with this plane. In both these approaches the universe is three 
dimensional and the background is expanding isotropically. 



The study of 2-D perturbations (like those due to pancakes, for example) in 
a 3-D expanding universe faces an operational problem: To begin with, we do 
not gain the dynamic range if we stick to 3D, even if we consider perturbations 
in a plane; the force between particles still has to computed by the solution of 
Poisson equation in three dimensions. Also, relevance of the interaction of matter 
outside the plane with these perturbations makes it, essentially, a 3D problem. 

Thus we are left with the second possibility. The two dimensional system is 
the intersection of an orthogonal plane and the "needles" and the force between 
the "particles" in this plane is given by the solution of the Poisson equation in 
two dimensions. Such a system is somewhat dichotomous with the background 
universe expanding isotropically. However, convenience is not the only reason for 
studying this somewhat strange system — relevant results for the evolution of 
density profiles around peaks in 2-D have also been computed for this type of a 
system [18]. 

Generalization of the NSR to the 2-D system was done using relations for 
cylindrical collapse by Padmanabhan [37] and we will test these predictions here. 

Although the system of infinite needles is appropriate for testing the pre- 
dictions in the intermediate regime, the same cannot be said for the asymptotic 
regime. We are dealing with a system that occupies a smaller number of dimen- 
sions in the phase space and the interaction of the constituents follows a different 
force law. Therefore, it is difficult to interpret, or carry over, results regarding 
stable clustering to the full 3-D system. 

4.1.1 Non-linear Scaling Relations 

The non-linear and the linear correlation functions at two different scales can be 
related by NSR. The relation between these scales is given by the characteristics 
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of the pair conservation equation [34]. For the two dimensional system of interest, 
this equation can be written as [39] 

dD dD 

— -h(A,x)— = 2h(A,X), (4.1) 

Here D = log(l+£), h = —v p /Hr is the scaled pair velocity, = 2x~ 2 f x r£(r)dr 
is the mean correlation function (£ is the correlation function), H is the Hubble's 
constant, X = log(x) and A = log(a). The characteristics of this equation are 
x 2 (l + £,(x, a)) = I 2 , where x and I are the two scales used in NSR. The self similar 
models due to Filmore and Goldreich [18] imply that for collapse of cylindrical 
perturbations the turn around radius and the initial density contrast inside that 
shell are related as x ta oc I /Si oc //£l(/). (Here £l is the linearly extrapolated 
mean correlation function). Noting that in two dimensions M oc x 2 , we find 
£(x) oc in the regime dominated by infall. Stable clustering limit implies 

£nl(o-,x) oc /) [39]. Thus in 2-D the scaling relations are 

^i(a, I) (Linear) 
£(a,x)oc< f L (a,/) 2 (Radial Infall) (4.2) 
£i(a, I) (Stable Clustering) 

A more general assumption compared to stable clustering involves taking h = con- 
stant asymptotically. In a system reaching steady state with both virialisation 
and mergers contributing to the evolution, one may reach a constant value for 
h, though it will not be unity if mergers are a dominant phenomenon. (This 
assumption has been discussed in, for example, [37].) It also allows a larger pa- 
rameter space to compare simulation results. If /i=constant asymptotically, then 
oc in this limit. Note that in 3D, the indices for three regimes are 1, 3 
and Sh/2 respectively. 

All features of clustering in three dimensions are present here as well. In 
particular, 



(i) If the asymptotic value of h scales with n such that h(n + 2) = constant 
then the final slope of the non-linear correlation function will be independent of 
the initial slope. 

(ii) If NSR exists then it will predict a specific index in the intermediate 
and asymptotic regimes which will depend on the initial power spectrum. In 
other words, existence of NSR implies that gravitational clustering does not erase 
memory of initial conditions. 

(iii) It is, however, possible that spectra which are not scale free acquire uni- 
versal critical indices at which the correlation functions grow in a 'shape invariant' 
manner. This comes about because the growth rate of correlation function varies 
with the local index and for an index that is not globally constant the correlation 
functions may 'straighten out' by this process. 

(iv) In 3-D clustering, n = — 1 in the intermediate regime and n = — 2 in the 
asymptotic regime [11] are the critical indices. These are the same for clustering 
in two dimensions. 

4.2 Simulations and Results 

We carried out a series of numerical experiments to test the ideas outlined above. 
We used a particle mesh code [9] to simulate power law models. The simulations 
were done with 1024 2 or 2048 2 particles in order to ensure that we had sufficient 
dynamic range to study all the three regimes in evolution of non-linear clustering. 
In particular, it is necessary to use larger simulations for power law spectra with 
a negative index. Here, we will present results for three models: n — 1, n — 
and n = —0.4. 

All the models are normalized by requiring the linearly extrapolated root 
mean square fluctuations in density, computed using a Gaussian filter, to be 
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unity at a scale of 10 grid points at a = 1.0. The results we present are for a — 1, 
2 and 5 for n = and n — 1, and a = 1, 2 and 3 for n = —0.4. 

A significant source of errors in large simulations is the addition of a small 
displacement in each step (fraction of a grid length) to a large position (up to 
2048 grid lengths). We avoid this problem by using net displacement for internal 
storage. 

We will show the correlation function and the pair velocity only for length 
scales larger than four grid lengths. We do this to avoid error due to shot noise 
and other artifacts introduced by various effects at smaller scales. This ensures 
that errors in our results are acceptably small. (Variations between different 
realizations give a dispersion of less than 10% in the correlation function.) 

In fig. 4. 2 we have plotted the non-linear correlation function £(x) as a function 
of the linearly extrapolated correlation function £l(0- Here the scales x and / 
are related by x 2 (l + £ ) = I 2 . Data for n — 1 is represented by circles, that for 
n = by stars and '+' marks the points for n = —0.4. Clearly, there are no 
systematic differences between the three models and the data points trace out 
a simple curve with three distinct slopes (We have also marked the 2 a errors 
calculated by averaging over several data sets. The error bars are plotted away 
from the NSR plot, for visibility and clarity). The NSR, shown as thick lines, is 



The slope in the intermediate regime is as expected. The asymptotic regime has 
a different slope than that predicted by stable clustering, which is shown as a 
dashed line. Unlike the observed relations for clustering in three dimensions, the 
coefficient for the intermediate regime is large. This has important implications 
for the critical index. 



Z(a,x) = l 2£ L (a,/) 2 
4.7&(M) 



3/4 



< 0.5; £(x) < 0.5 
0.5 < f £ (Z) < 2; 0.5 < £(x) < 8 

2<&(0; s<e» 



(4.3) 
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Figure 4.1: This figure shows the non-linear correlation function £(x) as a function 
of the linearly extrapolated correlation function Here the scales x and / are 

related by x 2 (l + £) = I 2 . Data for n — 1 is represented by circles, that for n = 
by stars and + marks the points for n = —0.4. For each of these models we have 
plotted data for the three epochs mentioned in the text. The estimated 2 a error bars 
are shown as vertical lines at three representative values of £ viz. at £=15.582, 1.65 
and 0.25, covering the nonlinear, intermediate and linear regimes. The error bars are 
shown away from the NSR plot for the sake of visibility. It is clear from this figure 
that there are no systematic differences between the three models and they trace out 
a simple curve with three distinct slopes. The slope of the curve in the intermediate 
regime is same as that predicted by the radial infall model. The stable clustering 
limit is shown as the dashed line and it is clear that the data points deviate from this 
curve. 



Panels of fig. 4.2 show as a function of x/x n i for the three models. These 
confirm that the slope of is consistent with the NSR shown in fig. 4.1. In 
each of these panels, the slope expected in the stable clustering limit is shown as 
a dashed line. 

As mentioned above, the existence of the NSR (eqn.(4.3)) implies that the 
slope of the correlation function will depend on the initial spectral index. To 
this extent, gravitational clustering does not erase memory of initial conditions. 
However, the differences of slope are significantly reduced by non-linear evolution. 



4.3 Summary 

Results obtained in this chapter can be summarized as follows: 

(i) We have verified that NSR for the correlation function exist for clustering 
in 2D in all the three regimes, just like in 3D. This NSR is independent of the 
power law index - at least for the three indices studied here. 

(ii) In the intermediate regime, the NSR in the form of eqn.(4.3), can be 
understood in terms of radial infall around peaks. Our simulations verify the 
predictions [39] for this regime. 

(iii) In the asymptotic regime, our results do not agree with the stable clus- 
tering hypothesis. The slope of the NSR in the asymptotic regime in fig. 4.1 
implies h = constant. We find that, in this regime, h ~ 3/4 for all the models 
studied here. 

(iv) The existence of NSR implies that the asymptotic slope of the correlation 
function depends on the initial slope. However, this is strictly true only for pure 
power law models; for other models it is possible for the spectra to be driven to 
a universal form. The NSR in the asymptotic regime seems to be linked to the 
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Figure 4.2: The correlation function as a function of rc/:r n ; for n = 1 model. 
Here rr n j oc a _2 ^ n+2 ^. Thick lines mark slopes expected from the non-linear scaling 
relations shown in fig. 4.1. The dashed line marks the expected slope of the cor- 
relation function in the stable clustering limit. The mismatch between the expected 
slope and the true slope in the intermediate regime may arise from the fact that the 
assumption of £ >> 1 used in computing the slope is not valid at the lower end of the 
regime. 




Figure 4.2: Similar plot for n — 0. 



logarithmic nature of the potential. 

A theoretical model for existence of such scaling relations, as mentioned be- 
fore, was derived on the basis of Spherical Collapse Model. In the next chapter 
we analyse the physical features of spherical collapse model in detail and derive a 
physically motivated approach to halting the collapse as opposed to the stardard 
'virialization' argument. 



Chapter 5 

An improved spherical collapse 
model 

Turbulence is life force. It is opportunity. Let's love turbulence 
and use it for change — Ramsay Clark 

5.1 Introduction 

Analytic modelling of the non-linear phase of gravitational clustering has been a 
challenging but interesting problem upon which a considerable amount of atten- 
tion has been bestowed in recent years. The simplest, yet remarkably successful, 
model for non-linear evolution is the Spherical Collapse Model (SCM, hereafter), 
which has been applied in the study of various empirical results in the gravita- 
tional instability paradigm. Unfortunately, this approach has serious flaws - 
both mathematically and conceptually. Mathematically, the SCM has a singular 
behaviour at finite time and predicts infinite density contrasts for all collapsed 
objects. Conceptually, it is not advisable to model the real universe as a sphere, 
in spite of the standard temptations to which theoreticians often succumb. The 
two issues are, of course, quite related, since, in any realistic situation, it is the 
deviations from spherical symmetry which lead to virialised stable structures get- 
ting formed. In conventional approaches, this is achieved by an ad hoc method 
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which involves halting the collapse at the virial radius by hand and mapping 
the resulting non-linear and linear overdensities to each other. This leads to the 
well known rule-of-thumb that, when the linear overdensity is about 1.68, bound 
structures with non-linear overdensities of about 178 would have formed. The 
singular behaviour, however, makes the actual trajectory of a spherical system 
quite useless after the turnaround phase — a price we pay for the arbitrary pro- 
cedure used in stabilizing the system. But the truly surprising feature is that, 
despite its inherent arbitrariness, the SCM, when properly interpreted, seems to 
give useful insights into the behaviour of real systems. The Press-Schechter for- 
malism [47], for the abundance of bound structures, uses SCM implicitly; more 
recently, it was shown that the basic physics behind the non-linear scaling rela- 
tions (NSR) obeyed by the two point correlation function can be obtained from 
a judicious application of SCM [37]. These successes, as well as the inherent 
simplicity of the underlying concepts, make the SCM an attractive paradigm for 
studying non-linear evolution in gravitational clustering and motivate one to ask: 
Can we improve the basic model in some manner so that the behaviour of the 
system after turnaround is 'more reasonable' ? 

It is clear from very general considerations that such an approach has to 
address fairly non-trivial technical issues. To begin with, exact modelling of 
deviations from spherical symmetry is quite impossible since it essentially requires 
solving the full BBGKY hierarchy. Secondly, the concept of a radius R(t) for a 
shell, evolving only due to the gravitational force of the matter inside, becomes 
ill-defined when deviations from spherical symmetry are introduced. Finally, our 
real interest is in modelling the statistical features of the density growth; whatever 
modifications we make to SCM should eventually tie up with known results for 
the evolution of, for instance, the two point correlation function. That is, we 



have to face the question of how best to obtain the statistical properties of the 
density field from the behaviour of a single system. 

In this chapter, these problems are addressed in a limited but focussed man- 
ner. The deviations from spherical symmetry are tackled by the retention of a 
term (which is usually neglected) in the equation describing the growth of the 
density contrast. Working in the fluid limit, we show that this term is physically 
motivated and present some arguments to derive an acceptable form for the same. 
The key new idea is to introduce a Taylor series expansion in (1/6) (where 8 is 
the density contrast) to model the non-linear evolution. We circumvent the ques- 
tion of defining the 'radius' of the non-spherical regions by working directly with 
density contrasts. Finally, we attempt to make the connection with statistical 
descriptors of non-linear growth, by using the non-linear scaling relations known 
from previous work. More precisely, we show that the modified equations predict 
a behaviour for the relative pair velocity (when interpreted statistically) which 
agrees with the results of N-body simulations. 

The chapter is divided into the following sections. The relevant equations 
describing the SCM are set out in Section 5.2; we also summarise the physical 
and ad hoc aspects of the SCM here. Next, we recast the equations in a different 
form and introduce two functions (i) a "virialization term" and (ii) a function 
hsc(8), whose asymptotic forms are easy to determine. The behaviour of h S c(8) 
in the presence and absence of the "virialization term" is also detailed here. In 
Section 5.4, we present the arguments that give the functional forms for the above 
term over a large range of 5; we then go on to present the results in terms of a 
single collapsing body and show how this term stabilizes a collapse which would 
have otherwise ended up in a singularity in terms of the growth of the density 
contrast with time. When this term is carried through into the equation for R(t) 



for a single system, it can be seen the radius reaches a maximum and gracefully 
decreases to a constant, remaining so thereafter. In the standard SCM, the 
radius decreases from the maximum all the way down to zero, thereby causing 
the density to diverge. Section 5.5 summarises the results and discusses their 
implications. 

5.2 The Spherical Collapse Model 

The scales of interest in the current work are much smaller than the Hubble 
length and the velocities in question are non-relativistic; Newtonian gravity can 
hence be used for the following analysis. We will consider the case of a dust- 
dominated, Q = 1 universe and treat the system in the fluid limit as being made 
up of pressureless dust of dark matter, with a smoothed density, p m (t,x), and a 
mean velocity, v(t, x). (This approach, of course, ignores effects arising from shell 
crossing and multi-streaming; these will be commented on later.) The density 
contrast, 5(x, t) is defined by 

p m (t,x)=p b (t)[l + 8(x,t)] (5.1) 

where p b denotes the smooth background density of matter. We define a velocity 
field u l = v 1 / (ad), where v 1 is the peculiar velocity (obtained after subtracting out 
the Hubble expansion) and a(t) denotes the scale factor. Taking the divergence 
of the field u % and writing it as 

diUj = Gij + e ijk Q k + ^5ij9 (5.2) 

where is the shear tensor, Q k is the rotation vector and 9 is the expansion, we 
can manipulate the fluid equations [38] to obtain the following equation for 5 



da 2 2a da 2a 2 



'd5\ 2 



+ (l + 5)(a 2 -2Q 2 ) (5.3) 



3(1 + 5) \da 

The same equation can be written in terms of time t as 
•■ 4 5 2 2d- 

AnG Pb 5(l + 5) + d 2 (l + 5) (a 2 - 2fi 2 ) (5.4) 

This equation turns out to be the same as the one for density contrast in the 
SCM, except for the additional term in (l + 5)(a 2 — 2f2 2 ), arising from the angular 
momentum and shear of the system. To see this explicitly, we introduce a function 
R(t) by the definition 

, . 9GMt 2 a 3 

l + ^=^-A- (5.5) 

where M and A are constants. Using this relation between S and R(t), equation 
(5.4) can be converted into the following equation for R(t) 

R=~- \a 2 (- 2 " 2^ 2 ) R (5.6) 

where the first term represents the gravitational attraction due to the mass inside 
a sphere of radius R and the second gives the effect of the shear and angular 
momentum. 

In the case of spherically symmetric evolution, the shear and angular momentum 
terms can be set to zero; this gives 

d 2 R GM 

IF = ~~kF (5 - 7) 

which governs the evolution of a spherical shell of radius R, collapsing under its 
own gravity; M can now be identified with the mass contained in the shell; this 
is standard SCM. 
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At this point, it is important to note a somewhat subtle aspect of these 
equations. The original fluid equations are clearly Eulerian in nature: i.e. the 
time derivatives give the temporal variation of the quantities at a fixed point in 
space. However, the time derivatives in equation (5.4), for the density contrast 
5, are of a different kind. Here, the observer is moving with the fluid element 
and hence, in this, Lagrangian case, the variation in density contrast seen by the 
observer has, along with the intrinsic time variation, a component which arises 
as a consequence of his being at different locations in space at different instants 
of time. When the 5 equation is converted into an equation for the function 
R(t), the Lagrangian picture is retained; in SCM, we can interpret R(t) as the 
radius of a spherical shell, co-moving with the observer. The mass M within 
each shell remains constant in the absence of shell crossing (which does not occur 
in the standard SCM for reasonable initial conditions) and the entire formalism 
is well defined. The physical identification of R is, however, not so clear in the 
case where the shear and rotation terms are retained, as these terms break the 
spherical symmetry of the system. We will nevertheless continue to think of R as 
the "effective shell radius" in this situation, defined by equation (5.5) governing 
its evolution. Of course, there is no such ambiguity in the mathematical definition 
of R in this formalism. 

Before proceeding further, let us briefly summarize the results of standard 
SCM. Equation (5.7) can be integrated to obtain R(t) in the parametric form 

r = ^i(i- cos e) (5.8) 

t = -^(9-sm9) (5.9) 

where Ri, 8^ and are the initial radius, initial density contrast and initial time, 
respectively, with Rf = (9GMi?/2)(l + 5^ ~ (9GMt 2 j2) for 6 t < I. Given 



M, there are only two independent constants, viz U and Si. All the physical 
features of the SCM can be easily derived from the above solution. Each spher- 
ical shell expands at a progressively slower rate against the self-gravity of the 
system, reaches a maximum radius and then collapses under its own gravity, 
with a steadily increasing density contrast. The maximum radius, R ma x = Ri/$i, 
achieved by the shell, occurs at a density contrast 5 = (97r 2 /16) — 1 ~ 4.6, which 
is in the "quasi-linear" regime. In the case of a perfectly spherical system, there 
exists no mechanism to halt the infall, which proceeds inexorably towards a sin- 
gularity, with all the mass of the system collapsing to a single point. Thus, the 
fate of the shell (as described by equations (5.8) and (5.9)) is to collapse to zero 
radius at 9 = 2n with an infinite density contrast; this is, of course, physically 
unacceptable. 

In real systems, however, the implicit assumptions that (i) matter is dis- 
tributed in spherical shells and (ii) the non-radial components of the velocities of 
the particles are small, will break down long before infinite densities are reached. 
Instead, we expect the collisionless dark matter to reach virial equilibrium. After 
virialization, \U\ = 2K, where U and K are, respectively, the potential and ki- 
netic energies; the virial radius can be easily computed to be half the maximum 
radius reached by the system. 

The virialization argument is clearly physically well-motivated for real sys- 
tems. However, as mentioned earlier, there exists no mechanism in the standard 
SCM to bring about this virialization; hence, one has to introduce by hand the as- 
sumption that, as the shell collapses and reaches a particular radius, say R max /2, 
the collapse is halted and the shell remains at this radius thereafter. This arbi- 
trary introduction of virialization is clearly one of the major drawbacks of the 
standard SCM and takes away its predictive power in the later stages of evo- 



lution. We shall now see how the retention of the angular momentum term in 
equation (5.6) can serve to stabilize the collapse of the system, thereby allowing 
us to model the evolution towards r vir = R max /2 smoothly. 

5.3 The hsc(8) function. 

As detailed in the previous section, the primary defect of the standard SCM is the 
ad hoc nature of the stabilization of the shell against its collapse under gravity, 
which arises on account of the assumption of perfect spherical symmetry, implicit 
in the neglect of the shear and angular momentum terms. We hence return to 
equation (5.3), retain the above terms, and recast the equation into a form more 
suitable for analysis. Using logarithmic variables, D sc = In (1 + 8) and a = In a, 
equation (5.3) can be written in the form (the subscript 'SC stands for 'Spherical 
Collapse') 

d 2 D sc 1 fdDscY | ldDsc 
da 2 3 \ da J 2 da 

^exp( J D S c)-l] + a 2 ( ( T 2 -2fi 2 ) (5.10) 

It is convenient to introduce the quantity, S, defined by 

S = a 2 (a 2 - 2fi 2 ) (5.11) 

which we shall hereafter call the "virialization term". The consequences of the 
retention of the virialization term are easy to describe qualitatively. We expect 
the evolution of an initially spherical shell to proceed along the lines of the stan- 
dard SCM in the initial stages, when any deviations from spherical symmetry, 
present in the initial conditions, are small. However, once the maximum radius 
is reached and the shell recollapses, these small deviations are amplified by a 
positive feedback mechanism. To understand this, we note that all particles in a 
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given spherical shell are equivalent due to the spherical symmetry of the system. 
This implies that the motion of any particle, in a specific shell, can be considered 
representative of the motion of the shell as a whole. Hence, the behaviour of the 
shell radius can be understood by an analysis of the motion of a single particle. 
The equation of motion of a particle in an expanding universe can be written as 

x, + 2- Xl = 1 5.12 

a a 1 

where a(t) is the expansion factor of the locally overdense "universe". The 
Xj term acts as a damping force when it is positive; i.e. while the background is 
expanding. However, when the overdense region reaches the point of maximum 
expansion and turns around, this term becomes negative, acting like a negative 
damping term, thereby amplifying any deviations from spherical symmetry which 
might have been initially present. Non-radial components of velocities build 
up, leading to a randomization of velocities which finally results in a virialised 
structure, with the mean relative velocity between any two particles balanced by 
the Hubble flow. It must be kept in mind, however, that the introduction of the 
virialization term changes the behaviour of the solution in a global sense and it is 
not strictly correct to say that this term starts to play a role only after recollapse, 
with the evolution proceeding along the lines of the standard SCM until then. It 
is nevertheless reasonable to expect that, at early times when the term is small, 
the system will evolve as standard SCM to reach a maximum radius, but will fall 
back smoothly to a constant size later on. 

The virialization term, S, is, in general, a function of a and x, especially 
since the derivatives in equation (5.4) are total time derivatives, which, for an 
expanding Universe, contain partial derivatives with respect to both x and t 
separately. Handling this equation exactly will take us back to the full non-linear 
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equations for the fluid and, of course, no progress can be made. Instead, we will 
make the ansatz that the virialization term depends on t and x only through 
<J(t,x). 

S(a,x) = S(6(a,x)) = S(Dsc) (5.13) 

In other words, S is a function of the density contrast alone. This ansatz seems 
well motivated because the density contrast, 8, can be used to characterize the 
SCM at any point in its evolution and one might expect the virialization term to 
be a function only of the system's state, at least to the lowest order. Further, the 
results obtained with this assumption appear to be sensible and may be treated 
as a test of the ansatz in its own framework. 

To proceed further systematically, we define a function /i sc by the relation 



dDsc 
da 



= 3/isc (5.14) 



For consistency, we shall assume the ansatz /isc( a ! x ) = ^sc [<H Q , X )]- The defi- 
nition of hsc allows us to write equation (5.10) as 

dh S C ,2 ^SC 1, / n x nl . S(D SC ) /rir\ 

= h sc- — +2 i e MDsc) - 1] + 3 (5.15) 

Dividing (5.15) by (5.14), we obtain the following equation for the function 

^sc(£>sc) 

= ^ - - + A- [exp(Dsc) - 1] + (5.16) 
dD sc 3 6 6/igc 9/isc 

If we know the form of either ^sc(-^sc) or S(D SC ), this equation allows us to 
determine the other. Then, using equation (5.14), one can determine -Dsc- Thus, 
our modification of the standard SCM essentially involves providing the form of 
S(Dsc) or hsc(Dsc)- We shall now discuss several features of such a modeling 
in order to arrive at a suitable form. 
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The behaviour of hsc(Dsc) can be qualitatively understood from our knowl- 
edge of the behaviour of 5 with time. In the linear regime (<J < 1), we know that 
S grows linearly with a; hence hsc increases with D$c- At the extreme non-linear 
end (5 1), the system "virializes" , i.e. the proper radius and the density of the 
system become constant. On the other hand, the density pb, of the background, 
falls like t~ 2 (or a~ 3 ) in a flat, dust-dominated Universe. The density contrast is 
defined by 5 = (p/pb) — 1 ~ p/ Pb (f° r <5 > 1) and hence 

5 oc t 2 oc a 3 (5.17) 

in the non-linear limit. Equation (5.14) then implies that hsc(8) tends to unity 
for (5 > 1. Thus, we expect that h S c(D sc ) will start with a value far less than 
unity, grow, reach a maximum a little greater than one and then smoothly fall 
back to unity. [A more general situation discussed in the literature corresponds 
to h — > constant as 5 — > oo, though the asymptotic value of h is not necessarily 
unity. Our discussion can be generalised to this case.] 

This behaviour of the hsc function can be given another useful interpretation 
whenever the density contrast has a monotonically decreasing relationship with 
the scale, x, with small x implying large S and vice-versa. Then, if we use a local 
power law approximation 5 oc x~ n for 5 ^> 1 with some n > 0, D sc oc ln(x _1 ) 
and 

gLDsc d Ml) xa v 
hsc oc — - — oc — — oc — oc -— (5.18) 
da a In a ax ax 

where v = ax denotes the mean relative velocity. Thus, hsc is proportional to 
the ratio of the peculiar velocity to the Hubble velocity. We know that this 
ratio is small in the linear regime (where the Hubble flow is dominant) and later 
increases, reaches a maximum and finally falls back to unity with the formation 
of a stable structure; this is another argument leading to the same qualitative 
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behaviour of the hsc function. 

Note that, in standard SCM (for which S = 0), equation (5.16) reduces to 

H dhsc ,2 ^sc 8 / r i r>\ 

3hsC dD^ c = ksc - ^ + 2 (5 ' 19) 

The presence of the linear term in 5 on the RHS of the above equation causes 
hsc to increase with 5, with hsc ^ S 1 ^ 2 for 5 > 1. If virialization is imposed as 
an ad hoc condition, then hsc should fall back to unity discontinuously — which 
is clearly unphysical; the form of S(S) must hence be chosen so as to ensure a 
smooth transition in hsc($) from one regime to another. 

As an aside, we would like to make some remarks on the nature of the "virial- 
ization term", S(S), in a somewhat wider context. As is well-known, gravitational 
clustering can be described at three different levels of approximation, by different 
mathematical techniques. The first approach tracks the clustering by following 
the true particle trajectories; this is what is done, for example, in N-body simulta- 
tions. This method does not involve any approximation (other than the validity 
of the Newtonian description at the scales of interest); it is, however, clearly 
analytically intractable. At the next level, one may describe the system by an 
one-particle distribution funtion and attempt to solve the collisionless Boltzmann 
equation for the distribution function f(t, x, v); the approximation here lies in 
the neglect of gravitational collisions, which seems quite reasonable as the time 
scale for such collisions is very large for standard dark matter particles. Finally, 
one can treat the system in the fluid limit described by five functions: the density 
p(t, x), mean velocity v(t, x), and gravitational potential <p(t, x), thus neglecting 
multi-streaming effects. Our analysis was based on this level of approximation. 
The key difference between the last two levels of description lies in the fact that 
the distribution function allows for the possibility of different particle velocities 
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at any point in space (i.e. the existence of velocity dispersions), while the fluid 
picture assumes a mean velocity at each point. It is also known that the gradi- 
ents in velocity dispersion can provide a kinetic pressure which will also provide 
support against gravitational collapse. While a detailed analysis of these terms 
is again exceedingly difficult, one can incorporate the lowest order effects of the 
gradient in the velocity dispersion by modifying equation (5.10) to the form 



where f(a,x) contains the lowest order contributions from the dispersion terms. 
We can then define 



and again invoke the ansatz S(a,x) = S(S). Note that S(S) now contains the 
lowest order contributions arising from shell crossing, multi-streaming, etc., be- 
sides the shear and angular momentum terms, i.e. it contains all effects leading 
to virialization of the system. It is explicitly demonstrated that velocity disper- 
sion terms arise naturally in the "force" equation (for the function h = —v/dx), 
derived from the BBGKY hierarchy, and play the same role as the function S(S) 
in the fluid picture. This clearly justifies the above procedure and shows that 
our approach could have a somewhat larger domain of validity than might be 
expected from an analysis based on the fluid picture. 




(5.20) 



S(a, x) = a 2 (a 2 - 2fi 2 ) + f(a, x) 



(5.21) 
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5.4 Detailed derivation of equation for h func- 
tion 

The zeroth and first moments of the 2 nd BBGKY equation [48] can be combined 
to obtain the following equation for the dimensionless function h = —v/ax 

3, i(1+?) a + ^-g- 9M »-- V = 
v ^' di 2 F A-kF 

din Fn dh\ j9 

h i( 4+ ^x-) + di- 2h - ^ 

where we have used the ansatz, h = h(£) ([19], [34]). In the above, we have 
defined 

£(x,a) = — J dx£(x,a)x 2 (5.23) 



F= dx +3{1 + i) ' x = lnx ( 5 - 24 ) 



1,1 ^■ l ' 2 7& ^ 



where II and £ are parallel and perpendicular peculiar velocity dispersions ([48], [27]). 
Finally, we have assumed that the 3-point correlation function has the hierarchi- 
cal form ([14], [48]) 

C123 = <2(£i 2 £i3 + 6s63 + 6263) (5.26) 

and defined 

/r 1 COS G 

dh[t(x) — (5.27) 

In the non-linear regime, £ 3> 1, the stable clustering ansatz yields a scale- 
invariant power-law behaviour for £ [14], with £ oc a^ 3_7 ^x~ 7 , if h — > 1 as £ — > 00. 
In this limit, we have 



^=(3-7)^ + 3 



(5.28) 
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and 



d\nF 

cW 



-7 



1 + 



3 -xir 



3- 7 /v r 



(5.29) 



2 

Further, we can write [57] M = M'x£ , where M' is a constant. Thus, equation 
(5.22) reduces, in the non-linear regime, to 

3 



, -dh , 9 h 3 

3 ^ -h 2 + -- 

d£ 2 3-7 



(4- 7 ) + 



37 



3MM- 

~4tT ]/ ~ 3-7 
dh 2 



+ 1 



hi 



OX ~ ^ + b 



(5.30) 



(3 - 7 )e. 

where we have retained terms upto order constant in £. We now assume that h 2 
and h\ are functions of £ alone, to first order. This yields 



, -dh , 9 h 3 

3 ^ - /i 2 + - - 

d£ 2 3-7 



3M ; r__ 
17 r 



3-7 



+ 1 



with 



c(0 



(4-7) + 



37 



(3 - 7)^ 

Clearly, if /i — > 1 as £ — > oo, we must have 



(5.31) 



(5.32) 



G(0 



3-7 



3M' r_ 



3 



3-7. 



+ 1 



'Mi 



(5.33) 



i.e. G(0 « -9M^/4tt(3 - 7) for £ > 1. 

Since G(£) tends to the above asymptote at late times, the residual part can be 
expanded in a Taylor series in l/£. Retaining the first two terms of the expansion 
in equation (5.31), we obtain 



-dh , 9 hl A B 3x 
3ft ^" ft + 2 + 2 = I + f + °« > 



(5.34) 



This is exactly the same as equation (5.40), with £ replacing 5. thus plays 

the same role as S{5) in the stabilising of the system against collapse. This 
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clearly implies that the velocity dispersion terms, h% and h\, will contribute to 
the support term; we are hence justified in writing the virialization term in the 
more general form 

S = a 2 (a 2 - 2fi 2 ) + f(a, x) (5.35) 

where f(a, x) contains contributions from effects arising from shell crossing, multi- 
streaming, etc. 

5.5 The virialization term 

We will now derive an approximate functional form for the virialization function 
from physically well-motivated arguments. If the virialization term is retained in 
equation (5.6), we have 

d 2 R GM H 2 R 

l^ = -l¥-— S (5 - 36) 
where H = a/a. Let us first consider the late time behaviour of the system. 
When virialization occurs, it seems reasonable to assume that R — ■> constant and 
R — > 0. This implies that, for large density contrasts, 

3r y M 

Using H = a/a = (2/3t), and equation (5.5) 

TJGMt 2 3, 3 



[n , • 7/ (<*»!) (5-38) 

Thus, the "virialization" term tends to a value of (—35/2) in the non-linear 
regime, when stable structures have formed. This asymptotic form for S(5) is, 
however, insufficient to model its behaviour over the larger range of density con- 
trast (especially the quasi-linear regime) which is of interest to us. Since S(5) 
tends to the above asymptotic form at late times, the residual part, i.e. the 
part that remains after the asymptotic value has been subtracted away, can be 
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expanded in a Taylor series in (1/5) without any loss of generality. Retaining the 
first two terms of expansion, we write the complete virialization term as 

S(S) = ~(1 + S) -j + ^ + OiS-") (5.39) 

Replacing for S(5) in equation (5.10), we obtain, for 5 ^> 1 

Q , r dhsc 2 h sc 1 A B 

3h5— - h sc + — + - = -j + - (5.40) 

[It can be easily demonstrated that the first order term in the Taylor series is 
alone insufficient to model the turnaround behaviour of the h function. We will 
hence include the next higher order term and use the form in equation (5.39) for 
the virialization term. The signs are chosen for future convenience, since it will 
turn out that both A and B are greater than zero.] In fact, for sufficiently large 5, 
the evolution depends only on the combination q = (B/A 2 ). This is most easily 
seen by rewriting equation (5.3), replacing S(5) with the above form. Taking the 
limit of large 5, i.e. 5 ^> 1, and rescaling 5 to 5 /A, we obtain 

d 2 5 3 d5 4 (d5\ 2 1 B 1 

db 2 + 2bdb ~ 37 \db) ~ ~^ 2 + A 2 ^ 2 5 ( ' ' 

= -~2+\ ( 5 - 42 ) 

From the form of the equation it is clear that the constants A and B occur in 
the combination q = B/A 2 and hence the non-linear regime is modelled by a one 
parameter family for the virialization term. 
Equation (5.36) can be written as 



•■ GM AR 
R = 



R 2 27t 2 



TJGMt 2 A B 

— V "ZTT 



(5.43) 



AR? 5 5 2 _ 

Using 5 = 9GMt 2 /2R 3 and B = qA 2 we may express equation (5.43) completely 
in terms of R and t. We now rescale R and t in the form R = r vir y(x) and 
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t = fix, where r vir is the final virialised radius [i.e. R — > r vir for t — > oo], where 
(3 2 = (8/3 5 )(j4/GM)rJ ir , to obtain the following equation for y(x) 



n V 4 27 y 7 

y =-i--?- 



(5.44) 



4 "x 6 

We can integrate this equation to find a form for y q (x) (where y q (x) is the function 
y(x) for a specific value of q) using the physically motivated boundary conditions 
y = 1 and = as x — > oo, which is simply an expression of the fact that 
the system reaches the virial radius r vir and remains here thereafter. The results 
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Figure 5.1: y q (x) for some values of q. The x axis has scaled time, x and the y axis 
is the scaled radius y. 

of numerical integration of this equation for a range of q values are shown in 
fig. 5.1. As expected on physical grounds, the function has a maximum and 
gracefully decreases to unity for large values of x [the behaviour of y(x) near 
x = is irrelevant since the original equation is valid only for 5 > 1, at least]. 
For a given value of q, it is possible to find the value x c at which the function 
reaches its maximum, as well as the ratio y max = Rmax/ r vir- The time, t max , at 
which the system will reach the maximum radius is related to x c by the relation 
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tmax = f3x c = t (l + z max )~ 3 / 2 , where t = 2/(3H ) is the present age of the 
universe and z max is the redshift at which the system turns around. Figure (5.2) 
shows the variation of x c and y max = (Rmax/r V i r ) for different values of q. The 
entire evolution of the system in the modified spherical collapse model (MSCM) 
can be expressed in terms of 

R(t) = r mr y q (t/(3) (5.45) 

where (3 = (t /x c )(l + z max )~ 3 / 2 . In SCM, the conventional value used for 
(r vir / R m ax) is (1/2), which is obtained by enforcing the virial condition that 
\U\ = 2K, where U is the gravitational potential energy and K is the kinetic 
energy. It must be kept in mind, however, that the ratio (r vir /R max ) is not really 
constrained to be precisely (1/2) since the actual value will depend on the final 
density profile and the precise definitions used for these radii. While we expect 
it to be around 0.5, some amount of variation, say between 0.25 and 0.75, cannot 
be ruled out theoretically. Figure (5.2) shows the parameter (R ma x/r V i r ), plotted 
as a function of q = B/A 2 (dashed line), obtained by numerical integration of 
equation (5.36) with the ansatz (5.39). The solid line gives the dependence of 
x c (or equivalently t max ) on the value of q. It can be seen that one can obtain 
a suitable value for the ij- vir / R max ) ratio by choosing a suitable value for q and 
vice versa. Using equation (5.14) and the definition 5 oc t 2 /R 3 , we obtain 

M*> = l-^| (5.46) 

which gives the form of hsc(%) for a given value of q; this, in turn, determines 
the function y q (x). Since 5 can be expressed in terms of x, y and x c as 5 = 
(97r 2 /2x 2 )x 2 /y 3 , this allows us to implicitly obtain a form for h S c(8), determined 
only by the value of q. Figure (5.3) shows the behaviour of hsc functions obtained 
by integrating equation (5.16) backwards, assuming that hsc ^ 1 as 5 ^ oo. It 
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Figure 5.2: The parameters (R ma x/ r vir) (broken line) and x c (solid line) as a function 
of q = B/A 2 . This clearly demonstrates that the single parameter description of the 
virialization term is constrained by the value that is chosen for the ratio r vir /R max . 
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Figure 5.3: The hsc function, obtained for various values of q. The values of q and 
y max = R m ax/ r vir f° r the curves are indicated at the top right hand corner. (Further 
discussion in text) 
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is seen that all the curves have the same turnaround behaviour expected on the 
basis of the physical arguments presented in the earlier section. 

If the functional form for hsc - determined, say, from N-body simulations - 
is used as a further constraint, we should be able to obtain the values of q. The 
major hurdle in attempting to do this is the fact that the available simulation 
results are given in terms of the averaged two point correlation function, £, and 
the averaged pair velocity, h(a,x), defined by 

£ = A ftix, a)x 2 dx ■ h(a, x) = - M^D (5.47) 
r 6 Jo ' ax 

where the two-point correlation function £ is defined as the Fourier transform 
of the power spectrum, P(k), of the distribution. The results published in the 
literature assume that h(a,x) depends on a and x only through £(a,x), that is, 
h(a,x) = h[£(a,x)]. This assumption has been invoked in several works in the 
past ([19], [34], [30], [37], [41]) and seems to be validated by numerical simula- 
tions. The fitting formula for h(£) can be obtained from related fitting formulas 
available in the literature [19]. These are, however, statistical quantities and are 
not well defined for an isolated overdense region. Hence we have to first make 
the correspondence between hsc (8) and h(£), which we do as follows. 
It is possible to show by standard arguments [34] that 

£ = 3h(l + £) (5.48) 



that is, 



^ = 3h (5.49) 
da 



where D — ln(l + £) and a — In a. Equation (5.49) is very similar to equation 
(5.14), which defines the function hsc(8), except for the different definitions of 
D and D sc in terms of £ and S respectively. This suggests that one can obtain a 
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relation between hsc(S) and by relating the density contrast 5 of an isolated 
spherical region to the two-point correlation function £ averaged over the distri- 
bution at the same scale. We essentially need to find a mapping between £ and 
5 which is valid in a statistical sense. 

Gravitational clustering is known to have three regimes in its growing phase, 
usually called "linear", "quasi-linear" and "non-linear" respectively. The three 
regimes may be characterized by values of density contrast as 5 1 in the linear 
regime, 1 < 5 < 100 in the quasi-linear regime and 100 < 5 in the non-linear 
regime. The three regimes have different rates of growth for various quantities 
of interest such as 5, £ and so on. In the linear regime, it is well known that 
the density contrast grows proportional to the scale factor, a. This implies that 
the power spectrum, P(k) = \Sk\ 2 (where 8k is the Fourier mode corresponding 
to 8(x)), grows as a 2 . Consequently, £, which is related to P(k) via a Fourier 
transform, also grows as a 2 , i.e. as the square of the density contrast. In the 
quasi-linear and non-linear regimes, the density contrast does not grow linearly 
with the scale factor and the relation between 8 and £ is not so clearly defined. 
The quasi-linear regime may be loosely construed as the interval of time during 
which the high peaks of the initial Gaussian random field have collapsed, although 
mergers of structures have not yet begun to play an important role. (This idea 
was used in [37] to model the non-linear scaling relations successfully). If we 
consider a length scale smaller than the size of the collapsed objects, the dominant 
contribution to £ (at this scale) arises from the density profiles centered on the 
collapsed peaks. Using the relation 




(5.50) 



for density profiles around high peaks, one can see that £ oc 8 in this regime. In 
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the non-linear regime, 5 and £ have the forms S(a,x) = a 3 F(ax), £ = a 3 G(ax), 
where x is a co-moving and r = arc is a proper coordinate. When the system 
is described by Lagrangian coordinates (which correspond to proper coordinates 
r = ax, i.e. at constant r), £ is proportional to 5. Thus, the relation £ oc 5 
appears to be satisfied in all regimes, except at the very linear end. Since we are 
only interested in the 5 > 1 range, we use £ ~ 5 and compare equations (5.48) 
and (5.14) to identify 

hsc(S)*h(0 (5.51) 

It is now straightforward to choose the value of q such that the known fitting 
function for the h function is reproduced as closely as possible. We use the 
original function given by Hamilton [19] to obtain the following expression for 

HO- 

m = 2 -(?* V ®Y (5.52) 

where V(£) is given by the fitting function 

= 7 ( 1 + 0.0158 £ 2 + 0.000115 £ 3 \ 1/3 
lU ? ^1 + 0.926 £ 2 - 0.0743 £ 3 + 0.0156 £ 4 ; l ' ' 

Figure (5.4) shows the simulation data represented by the fit (solid line) [19] 
and the best fit (dashed line), obtained in our model, for q ~ 0.02. We note that 
the fit is better than 5% for all values of density contrast 5 > 15. The change 
in the fit is very marginal if one imposes the boundary condition h(5) — > 1 for 
5 ^> 1, instead of constraining the curves to match at their peaks (for example, 
the change in the peak height is ~ 1%, if we impose the above condition at 
5 = 10000). 

Figure (5.5) shows the plot of scaled radius y q (x) vs x, obtained by integrating 
equation(5.44), with q = 0.02. The figure also shows an accurate fit (dashed line) 
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Figure 5.4: The best fit curve for the h function (dashed line) to the simulation data 
(solid line). The simulation results are obtained from Hamilton [19] and the fit is 
obtained by adjusting the value of q parameter until the curves coincide. 
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Figure 5.5: Plot of the scaled radius of the shell y q as a function of scaled time 
x (solid line) and the fitting formula y q — (x + ax 3 + bx 5 )/(l + cx 3 + bx 5 ), with 
a = —3.6, b = 53 and c = —12 (dashed line) (See text for discussion) 
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to this solution of the form 



Vq( X ) 



x + ax 3 + bx z 



(5.54) 



1 + cx 3 + bx 5 

with a = —3.6, b = 53 and c = —12. This fit, along with values for r V i r and 
Zmax, completely specifies our model through equation (5.45). It can be observed 
that (r V i r /R max ) is approximately 0.65. It is interesting to note that the value 
obtained for the [r vir j R max ) ratio is not very widely off the usual value of 0.5 used 
in the standard spherical collapse model, in spite of the fact that no constraint 
was imposed on this value, ab initio, in arriving at this result. Part of this 
deviation may also originate in the fit which has been used for Hamilton et 
al. [19] noticed that objects virialised at R ma x/ r vir ~ 1-8, instead of 2, in their 
simulations. 
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Figure 5.6: The non-linear density contrast in the SCM (solid line) and in the modified 
SCM (dashed line), plotted against the linearly extrapolated density contrast 5l- 



Finally, figure (5.6) compares the non-linear density contrast in the modified 
SCM (dashed line) with that in the standard SCM (solid line), by plotting both 
against the linearly extrapolated density contrast, 5^. It can be seen (for a given 



system with the same z max and r vir ) that, at the epoch where the standard SCM 
model has a singular behaviour (Sl ~ 1.686), our model has a smooth behaviour 
with 5 ~ 110 (the value is not very sensitive to the exact value of q). This is not 
widely off from the value usually obtained from the ad hoc procedure applied in 
the standard spherical collapse model. In a way, this explains the unreasonable 
effectiveness of standard SCM in the study of non-linear clustering. 

As mentioned earlier, deviations from spherical symmetry are expected to be 
small at early epochs and to grow as the system evolves. One would thus expect 
the two curves of figure (5.6) to approach each other as S — > 1 (from above). 
Further, the curves should overlay in the linear regime (Sl < 1). It can be seen 
from the figure that the 2 curves do approach each other as Sl reduces towards 
unity. However, the MSCM has been obtained using a Taylor expansion in (l/S); 
it is clearly not applicable for 5 <C 1. Further, the region S > 15 has been used 
to fit the function h(S) to the data of Hamilton [19]. Hence, one cannot compare 
the curves in the linear regime. 

Figure (5.6) also shows a comparison between the standard SCM and the 
MSCM in terms of S values in the MSCM at two important epochs, indicated 
by vertical arrows, (i) When R = R max /2 in the SCM, i.e. the epoch at which 
the SCM virializes, 5(MSCM) ~ 83. (ii) When the SCM hits the singularity, 
(S L ~ 1.6865), 5(MSCM) ~ 110. 

We note, finally, that figure (5.6), which shows the effects of evolution as a 
mapping from linear to non-linear density contrasts, contains a subtle implicit 
assumption regarding the definition of the non- linear density contrast. The radius 
R of a system is not a rigorously defined quantity in the absence of spherical 
symmetry, and obviously, any argument involving 'virialization' precludes strict 
spherical symmetry. It is, however, a conventional practice to define the 'radius', 



R, even for a virialized system without strict spherical symmetry. For example, 
this approach is used to define the density contrast at 'virialization' (which has 
the value, S V i r ~ 200) in the standard SCM. In our model we have an explicit 
equation for R; once R and M are given, the non-linear density contrast is a 
well-defined quantity. 

5.6 Results and Summary 

It has been shown how the Taylor expansion of a term in the equation for the 
evolution of the density contrast, 5, in inverse powers of 5, allows us to have a 
more realistic picture of spherical collapse, which is free from arbitrary abrupt 
"virialization" arguments. Beginning from a well motivated ansatz for the de- 
pendence of the "virialization" term on the density contrast we have shown that 
a spherical collapse model will gracefully turn around and collapse to a constant 
radius with 5 ~ 110 at the same epoch when the standard model reaches a sin- 
gularity. Figure (5.5) shows clearly that the singularity is avoided in our model 
due to the enhancement of deviations from spherical symmetry, and consequent 
generation of strong non-radial motions. 

We derive an approximate functional form for the virialization term starting 
from the physically reasonable assumption that the system reaches a constant 
radius. This assumption allows us to derive an asymptotic form for the virializa- 
tion term, with the residual part adequately expressed by keeping only the first 
and second order terms in a Taylor series in (1/5). It is shown that there exists 
a scaling relation between the coefficients of the first and second order terms, 
essentially reducing the virialization term to a one parameter family of models. 

The form of the h function published in the literature, along with a tentative 
mapping from S to £, in the non-linear and quasi-linear regimes, allow us to further 



constrain our model, bringing it in concordance with the available simulation 
results. Further, it is shown that this form for the virialization term is sufficient 
to model the turnaround behaviour of the spherical shell and leads to a reasonable 
numerical value for density contrast at collapse. 

There are several new avenues suggested by this work. 

(i) The assumption h$c — > 1, R — > r vir is equivalent to "stable clustering", 
in terms of the statistical behaviour. Since stable clustering has so far not been 
proved conclusively in simulations and is often questioned, it would be interesting 
to see the effect of changing this constraint to h — > constant for t — > oo. 

(ii) The technique of Taylor series expansion in (1/5) seems to hold promise. It 
would be interesting to try such an attempt with the original fluid equations and 
(possibly) with more general descriptions. 

(hi) It must be stressed that we used the 5 — £ mapping — possibly the weakest 
part of our analysis, conceptually — only to fix a value of q. We could have 
used some high resolution simulations to actually study the evolution of a real- 
istic overdense region. We conjecture that such an analysis will give results in 
conformity with those obtained here. 

(iv) Finally, the curves of figure (5.6) can be used to describe the spatial distri- 
bution of virialised haloes ([30], [52]). It would be interesting to investigate how 
things change when the MSCM is used in place of the standard spherical collapse 
model. 

In the next chapter we shall approach the problem of structure formation 
from a completely different point of view, namely that of Quasi Steady State 
Cosmology (QSSC). 



Chapter 6 

Structure formation in QSSC 



Just because everything is different doesnt mean anything 
has changed - Paradoxial quotations 

6.1 Introduction 

The Quasi-Steady State Cosmology (QSSC) was first proposed in 1993 and ex- 
plored further by Fred Hoyle, Geoffrey Burbidge and Jayant Narlikar in a series 
of papers ([20], [21], [22], [23], [24]). The QSSC offers an alternative to the com- 
monly accepted big bang cosmology, and the above work claims to provide a 
singularity-free cosmological model, which is consistent with the data on discrete 
source populations, can explain the production of light nuclei as well as the spec- 
trum and anisotropy of the microwave background. Because the dynamical and 
physical conditions in this cosmology are considerably different from those in the 
standard cosmology, the theoretical reasoning required to understand what is 
observed may differ too. In short, one may not simply lift a theoretical line of 
reasoning from standard cosmology and expect to apply it to the same problem 
in the QSSC. 

One of the outstanding problems in modern big bang cosmology is the prob- 
lem of formation of large scale structure in the universe. The standard approach 
consists in starting with prescribed primordial fluctuations of spacetime geome- 
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try and matter density, evolving them through an inflationary era, having them 
interact with nonbaryonic dark matter, then carrying out N-body simulations of 
interacting masses which may eventually form into groups to be identified with 
large scale structures like galaxies, clusters, superclusters and voids, etc. Al- 
though a lot of this work has gone into cosmology textbooks ([45], [36]), it is a 
fair comment to say that no unique and generally acceptable structure formation 
scenario has yet emerged in standard cosmology. 

The problem of structure formation poses a challenge in the QSSC also and 
it should be viewed against the background of the above standard approach. As 
we shall see in Section 6.2, the QSSC does not have an era when the baryonic 
matter density in the universe was ~ 10 81 times its present value, as it was in 
the big bang cosmology in the immediate post-inflation era. Thus the growth of 
fluctuations in the form of gravitational instabilities will not be similar in this 
cosmology to that in the big bang cosmology. 

Recently the gravitational stability of the QSSC models against small pertur- 
bations was examined in detail in a paper by Banerjee and Narlikar [12]. They 
found the cosmological solution to be stable and thus there was no net growth 
in density fluctuations. The model is basically oscillatory and perturbations of 
density and metric grow only to a finite amount during the contraction phase 
and then decay during the expansion phase. These authors concluded that grav- 
itational instability alone cannot lead to formation of structures in the QSSC. 
Instead, explosive matter creation in the so-called minibangs is expected to be 
the principal cause of forming structures. In this chapter we try to understand 
the pattern of formation and growth of structures in the QSSC though numerical 
simulations by using a simplified toy model. 

The organization of this chapter is as follows: In section 6.2 we briefly review 



the basic theory of QSSC. The numerical toy model will be introduced in section 
6.3. Section 6.4 is devoted to computing the two point correlation function for the 
distributions arising in the toy model and its comparison with observations. In 
section 6.5 we summarise the results by highlighting the success of this approach 
and indicating how it can be further improved. 

6.2 The Basic Theory of the QSSC 

The basic formulation of the QSSC is via the Machian theory of gravity first 
proposed by Hoyle and Narlikar ([25], [26]) in which the origin of inertia is linked 
to a long range scalar interaction between matter and matter. Specifically, the 
theory is derivable from an action principle with the simple action: 

A = ~Y, / m ads a , (6.1) 

where the summation is over all the particles in the universe, labeled by the index 
a, the mass of the ath particle being m a . The integral is over the world line of 
the particle, ds a representing the element of proper time of the ath particle. 

The mass itself arises from interaction with other particles. Thus the mass 
of particle a at point A on its worldline arises from all other particles b in the 
universe: 

m a = J2 m (b)(A), (6.2) 

where m^){X) is the contribution of inertial mass from particle b to any particle 
situated at a general spacetime point X. The long range effect is Machian in 
nature and is communicated by the scalar mass function m^) (X) which satisfies 
the conformally invariant wave equation 

Dm (b) + \ Rm (b) + m \b) = N (b)- (6-3) 
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Here the wave operator is with respect to the general spacetime point X. R 
is the scalar curvature of spacetime and the right hand side gives the number 
density of particle b. The field equations are obtained by varying the action 
with respect to the spacetime metric g^. The important point to note is that 
the above formalism is conformally invariant. In particular, one can choose a 
conformal frame in which the particle masses are constant. If the constant mass 
is denoted by m p , the field equations reduce to 

Rik _ }_ g ik R + Ag ik = _^f_{ T ik _ f( C i C k _ ]_gik c i Ci ^ ( 6 4 ) 

where c is the speed of light and C is a scalar field which arises explicitly from 
the ends of broken world lines, that is when there is creation (or, annihilation) 
of particles in the universe. The constant / denotes the coupling of the C-field 
to spacetime. Thus the divergence of the matter tensor T lk need not always be 
zero, as the creation or annihilation of particles is compensated by the non-zero 
divergence of the C-field tensor in Eq.(6.4). The quantities G (the gravitational 
constant) and A (the cosmological constant) are related to the large scale distri- 
bution of particles in the universe. Thus, 

G = 4^f' A = ~iv%5' (6 - 5) 

N being the number of particles within the cosmic horizon. 

Note that the signs of the various constants are determined by the theory 
and not put in by hand. For example, the constant of gravitation is positive, the 
cosmological constant negative and the coupling of the C-field energy tensor to 
spacetime is negative. 



6.2.1 Matter Creation 

The action principle tells us that matter creation is possible at a given spacetime 
point provided the ambient C-field satisfies the equality CjC* = w? p at that point. 
In normal circumstances, the background level of the C-field will be below this 
level. However, in the strong gravity obtaining in the neighbourhood of compact 
massive objects, the value of the field can be locally raised. This leads to creation 
of matter along with the creation of negative C-field energy. The latter also has 
negative stresses which have the effect of blowing the spacetime outwards (as in 
an inflationary model) with the result that the created matter is thrown out in an 
explosion. Qualitatively, the creation and ejection proceeds along the following 
lines. 

The process normally begins by the creation of the C-field along with matter 
in the neighbourhood of a compact massive object. The former, being propagated 
by the wave equation, tends to travel outwards with the speed of light, leaving 
the created mass behind. However, as the created mass grows, its gravitational 
redshift begins to assert itself, and the C-field gets trapped in the vicinity of the 
object. As its strength grows, its repulsive effect begins to manifest itself, thus 
making the object less and less bound and unstable. Finally, a stage may come 
when a part of the object is ejected from it with tremendous energy. It is thus 
possible for a parent compact mass to eject a bound unit outwards. This unit 
may act as a center of creation in its own right. 

We shall refer to such pockets of creation as minibangs or mini- creation events 
(MCEs). A spherical (Schwarzschild type) compact matter distribution will lead 
to a spherically symmetric explosion whereas an axi-symmetric (Kerr type) dis- 
tribution would lead to jet-like ejection along the symmetric axis. Because of the 
conservation of angular momentum of a collapsing object, it is expected that the 
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latter situation will in general be more likely. 

In either case, however, the minibang is nonsingular. There is no state of 
infinite curvature and terminating worldlines, as in the standard big bang, nor is 
there a black hole type horizon. The latter because the presence of the C-field 
causes the collapsing object to bounce outside the event horizon. 

6.2.2 The Cosmological Solution 

The feedback of such minibangs on the spacetime as a whole is to make it expand. 
In a completely steady situation, the spacetime will be that given by the de Sitter 
metric. However, the creation activity passes through epochs of ups and downs 
with the result that the spacetime also shows an oscillation about the long term 
steady state. Sachs [49] has computed the general solutions of this kind and the 
simplest such solution with the line element given by 

ds 2 = c 2 dt 2 - S 2 (t)[dr 2 + r 2 (d6 2 + stn 2 6d(j) 2 )}, (6.6) 

where c stands for the speed of light has the scale factor given by 

27rr(t)" 



S(t) = e l ' p 



1 + rj cos 



(6.7) 



Q 

The constants P and Q are related to the constants in the field equations, while 
r(t) is a function ~ t which is also determined by the field equations. We shall, 
however, use the approximation r(t) = t which is adequate for the approach used 
in this chapter. The parameter rj may be taken positive and is less than unity. 
Thus the scale factor never becomes zero: the cosmological solution is without a 
spacetime singularity. The form of the scale factor, S(t), in the metric (6.6) is 
shown in Figure 6.1. 




(t/Q) 



Figure 6.1: The scale factor S(t) of the QSSC in the upper panel against t to show 
how several oscillatory cycles of short period Q are accommodated in the longer e- 
folding time P of the exponential expansion. In the lower panel are sketched a few 
oscillations on an expanded timescale with our present epoch marked. 

6.2.3 Observational Checks 

Hoyle ([21], [22]) have shown that the above cosmology gives a reasonably good fit 
to the observations of discrete source populations, such as the redshift-magnitude 
relation, radio source count, angular diameter- redshift relation and the maximum 
redshifts so far observed, with the choice of the following set of parameters: 

P^20Q, Q^AA x 10 10 yrs, r] = 0.8, A = -0.3 x 10~ 56 cm" 2 , t = 0.7Q. 

Of these, the last is the present epoch of observation. It is not essential that 
the model should have only these parameteric values. Indeed, the parameter space 
is wide enough to make the model robust. Moreover, the fitting of observations to 
theory does not require postulating ad hoc evolution which is commonly necessary 
in the case of standard cosmology. 

The above framework thus outlines a cosmological model without a beginning 
and without an end, in which a de Sitter type exponential expansion, character- 



ized by a very long time scale P, is superposed with finite size oscillations of a 
shorter time scale Q. Each cycles are statistically identical in their physical prop- 
erties. In this sense the universe is 'quasi-steady'. We next see how structures 
might grow and proliferate in such a universe. 

6.3 A Toy Model for Formation of Structures 

In an attempt to understand how structures may possibly grow and distribute in 
space we have carried out the following numerical experiment in two as well as 
three dimensional space. We describe the 2D case first and detail the 3D versions 
subsequently. 

6.3.1 2D-Simulations 

A large number of points (N ~ 10 5 — 10 6 ), each one representing a mini-creation 
event, is distributed randomly over a unit square area . The average nearest- 
neighbour distance for such a distribution will then be (1/y/N). Now suppose 
that in a typical mini-creation event, each particle generates another neighbour 
particle at random within a distance, d = x/y/N in 2D . Here, the number x 
is a fraction between and 1 . We shall call x the separation parameter. As 
explained in section 6.2.1, the above denotes an ejected piece lying at a distance 
< d from the original compact object. 

The sample area is then uniformly stretched by a linear factor \/2 to represent 
expansion of space. We now have the same density of points as before, i.e., 2N 
points over area of 2 units. From this enlarged square remove the periphery so 
as to retain only the inner unit square. This process thus brings us back to the 
original state but with a different distribution of an average N points over a unit 
square. This process is repeated n times . Here the number of iterations, n, plays 
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the role of "time" as in the standard models of structure formation. The number 
distribution of points evolves as the 'creation process' generates new points near 
the existing ones. We will refer to each point as a 'particle' or 'unit'. 

Not surprisingly, soon after, i.e., after n = 3 — 4 iterations of the above 
procedure, clusters and voids begin to emerge in the sample area and create 
a Persian Carpet type of patterns. As the experiment is repeated, voids grow 
in size while clusters become denser. Figure 6.2 illustrates a typical numerical 
simulation. It shows that expansion coupled with creation of matter is a natural 
means of generating voids and clusters. But what of the filaments ? Here we 



Figure 6.2: A cluster-void distribution generated in the 2D toy model for iV = 
100,000 initially randomly distributed particles, with typical separation parameter 
x = 0.8, and the number of iterations n = 10. Each particle resembles a galaxy. For 
further discussion see the text. 

recall that the creation process near a typical compact massive object will not be 
isotropic if the mass is spinning. Matter will be preferentially ejected along the 
axis of spin. To build this effect into the above simulation we adopt the following 
algorithm. 



N = 100,000, x = 0.8, n = 10 




o.o 

X DIRECTION 



We assume that in a typical n > 2 iteration, the creation of the new neighbour 
unit C around a typical unit B is not entirely random, but, instead, related to 
the previous history of creation of B from an earlier generation unit A. So the 
direction BC is broadly aligned with the direction AB in which B itself was 
ejected. Typically this is ensured by assuming that the ejection is at a random 
angle in the forward semicircle as explained in Figure 6.3. We will refer to this 
as aligned ejection, as opposed to the isotropic ejection of Figure 6.2. Physically 




Figure 6.3: Schematic procedure for creating units in aligned direction for n > 2. 
Particle A is a representative of first generation units which are distributed randomly. 
B is a representative of the second generation units, being created in a random 
direction. Point C represents a third generation unit which has been created in the 
half plane lying away from A off the line perpendicular to AB. BC therefore makes 
an acute angle, ip, with the line ABB' . 

this means that the unit B ejected by A retains 'memory' of its origin through 
its spin which is more or less aligned with the spin of A. Which is why when it 
ejects a unit C, it is more or less aligned with the earlier ejection direction AB. 

Although this algorithm does not demand strict alignment, it is interesting 
to note that the filamentary structure grows along with voids as n increases. 
Features generated in this way have very suggestive similarities with the observed 
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large scale structure as shown in a typical simulation of Figure 6.4. We have 
also investigated the result of restricting the secondary ejection to a narrower 
angle, e.g, by keeping the angle f of Figure 6.3 in the range (— 7r/4, 7r/4). Not 
surprisingly we find the filamentary structure more pronounced in such a case. 
In general, we may argue that the higher the angular momentum per unit mass 
of the compact object causing ejection, the narrower is the angle of ejection, the 
greater is the alignment and hence more pronounced the filamentary structure. 
Figure 6.5 shows a typical two dimensional gravitational clustering simulation 

N = 100,000, x = 0.8, n = 10 

0.2 

Z 
o 

H 

-0.2 



-0.4 



-0.4 -0.2 0.0 0.2 0.4 

X DIRECTION 

Figure 6.4: A computer simulated filament-void distribution with n(> 2) iterations 
of aligned ejections having new points following the rule of Figure 6.3, for the same 
parameters of Figure 6.2. 

data in standard big bang cosmology. 

It can be seen that both compact and extended structures are present in both 
the approaches to structure formation. Since there is no observational data with 
which comparisons can be made in two dimensions we shall henceforth deal with 
the 3D simulations only. 
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Figure 6.5: A power law 2D simulation in the standard big bang cosmology for power 
index, n = —0.4 of density fluctuations. 

6.3.2 3L>-Simulations 

The 3D simulation is similar, with the necessary modifications for the higher 
dimensionality. Thus we start with a unit cube with iV points distributed at 
random within it, the typical interpoint distance being (1/ \/iV). Creating a new 
near neighbour for each particle by the same rule as in the 2D case, we need to 
expand each edge of the cube by the factor \/2. We next apply the same algo- 
rithm favouring aligned ejection, suitably modified for 3D. To compare the three 
dimensional distributions with the observed distributions made up from redshift 
surveys, we need to take a thin inside slice of the cube perpendicular to one of 
its edges and examine the distribution of points therein. Before making such a 
comparison, however, we will first apply 3-D-simulations within the framework of 
the QSSC. 



6.3.3 Simulations of QSSC Cycles 

To bring the toy model closer to the reality of the QSSC, we proceed as follows. 
We expect the creation activity to be confined largely to a narrow era around a 
typical oscillatory minimum, when the C-field is at its strongest. By considering 
the number density of collapsed massive objects at one oscillatory minimum of 
QSSC to be /, the number density at the next oscillatory minimum would fall to 
f exp (-3Q/P), if no new massive objects were added. Thus to restore a steady 
state from one cycle to the next, 

af=[l- exp (-3Q/P)]/ ~ (3Q/P)/, (6.8) 

masses must be created anew. In other words, a fraction (3Q/P) of the total 
number of massive objects must duplicate themselves in the above fashion. 

Notice that, unlike the old steady state theory which had new matter appear- 
ing continuously, we have here discrete creation, confined to epochs of minimum 
of scale factor. The 'steady-state' is maintained from one cycle to next. Which is 
why the above addition af is required at the beginning of each cycle. Therefore, 
instead of creating a new neighbour particle around each and every one of the 
original set of N particles, we do so only around aN of these points chosen ran- 
domly, where the fraction a is as defined in (6.8). Likewise, the sample volume 
is homologously expanded by the factor exp (3Q/P) only instead of by factor 2. 
We choose the inner cube as before. Figure 6.6 shows the simulated distributions 
in cubical slices for isotropic as well as aligned ejections. After a few iterations 
clusters and voids begin to appear, with the case for aligned ejections showing 
filaments. For a comparison, see an actually observed distribution of galaxies 
from a redshift survey in Figure 6.7. In the above approximation we have as- 
sumed that the creation activity is concentrated at the oscillatory minima. It 
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Figure 6.6: A 3-dimensional version, adapted for the QSSC with iV = 1,000,000, 
x = 0.3, n = 10, and P/Q = 20. Slice thickness in the Z direction is AZ = 0.001. 
Evidently, voids are seen separated by filamentary structures. The left panel is for the 
case of isotropic distribution of particles, whereas the right panel shows the case of 
aligned ejections. 

could be extended over an appreciable part of the oscillatory period, in which 
case one would see large scale structure in the radial direction as seen from an 
observer. We have not modified an algorithm to cover such cases, but feel that 
this should be investigated, especially since the recent analysis of the redshift- 
magnitude relation for supernovae has generated interest in the QSSC models of 
this kind [13]. 

6.4 The Two Point Correlation Functions 

Although visual inspection of Figures 6.6 and 6.7 suggests that the simulation 
is proceeding along the right lines, a quantitative measure of the cluster-void 
distribution will help in comparing simulations with reality. 
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Figure 6.7: A FLAIR redshift survey in the direction of Hydra-Centaurus[54] 
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The dimensionless autocorrelation function 



£(r) =< [p(r)- < p >][p(n + r)- < p >] > / < p > 2 , 



(6.9) 



where < p > is the average density in the volume, is one convenient measure of 
such irregularities in the space distribution. Typically, different classes of objects 
cluster at different characteristic lengths. To fix ideas in the present model, 
we will look at distribution of clusters of galaxies. Observationally, it is believed 
that the two point correlation function for cluster distribution obeys the following 
scaling law: 



with 7 ~ 1.8 and r = 25k" 1 Mpc, where the Hubble constant at the present 
epoch is taken to be lOO/i kms _1 Mpc _1 . In order to quantify the issues of forma- 
tion of structures in this scenario we have taken the following measures. 

It is known that instead of having a smooth distribution of matter on large 
scales, the observed universe has structures of typical sizes of a few tens of mega- 
parsecs . These "structures" are regions of density considerably higher than the 
background density, with the maximum density contrast 5 = (5p(r)/ < p >) 
going from order unity (in the case of clusters) to a few thousand (in the case of 
the galaxies). 

Any process which generates structures must be able to produce to the ze- 
roth order , entities whose density contrast is of such magnitude and with the 
property that on larger and larger distance scales, the density contrast becomes 
less significant. This is to ensure that on a large enough scale the universe is 
homogeneous. 

Given this prescription for generating structures without gravitational dy- 
namics, we first ensured that the visual impression created by the cluster sim- 




(6.10) 



ulations did imply that as the number of iterations were increased the number 
of high density regions also increased. In the initial configuration one expects to 
find regions of high density arising only because of the Poisson noise. In the later 
"epochs" after a few iterations, however, we expected and did find that the varia- 
tion of the one point distribution function for density (p/ < p >) with < p > the 
average density in the volume, showed a steady and significant increase in the the 
number of high and intermediate density regions, as is expected in a clustering 
scenario. We further observed that the value of maximum density also increased 
as a function of the number of iterations, which in this experiment corresponds to 
"time" . The density field has been generated on a grid placed into the simulation 
volume using the algorithm of cloud in cell. 

Our simulations show the growth of structures through rise in the density 
maximum as a function of number of iterations. The aligned ejection mode leads 
to faster clustering than the isotropic one. 

One must also examine the dependence of this "growth" on another important 
parameter in this prescription, namely the typical maximum separation between 
a creation site and the unit which is created. This was indicated by the parameter 
x in our earlier discussion. 

Again our studies investigate results of the structure formation algorithm 
when the parameter x is changed. We find that higher densities are achieved 
when this distance (in units of boxsize) is made smaller. This is intuitively 
expected. 

In the QSSC case, clustering is stronger in the early epochs for the isotropic 
ejection model, although at a later stage the density function for the aligned 
ejection model catches up and ultimately exceeds the rate for the isotropic case. 

The next quantitative measure that we computed from these data set was 
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the two point correlation function. The following figures summarise the results 
of these computations. It can be seen that the observationally obtained power law 
dependence of the two point correlation function £(r) = (r /r) L8 can be obtained 
provided a sufficient number of iterations has been performed, i.e, a sufficient 
length of "time" has elapsed. 

Figure 6.8 shows the two point correlation function for the case of the QSSC 
based model. As "time" goes on, the slope of the correlation function gets closer 
and closer to —1.8. From the value of the X-axis intercept of the two point corre- 
lation function it is possible to get a rough estimate of the size of the structures 
in units of the size of our simulation box. From our results we have estimated 
that the size of the structures formed is approximately (3 = 0.15 — 0.3 times the 
boxsize. If one sets these values equal to the observationally accepted value of 
r , one can get a better physical sense of the results. If we set, f3 = 0.3, say, and 
r = 25/iT 1 Mpc, then the linear size of the simulation box would be ~ 84/i _1 
Mpc. 

The above exercise is an attempt to relate our toy model to a realistic cos- 
mological scenario. The model per se talks of a 'dimensionless' box containing N 
points. With the above identification, we have 10 5 points in a volume of (84/iT 1 ) 3 
(Mpc) 3 . Let us assign a mass of 10 n M Q to each point. We then get a cosmological 
smoothed out density of 

p = , ; , " , xq = 0.6/i x 10 n - 12 Pc (6.11) 
1 (84/i- 1 ) 3 (Mpc) 3 H y J 

where p c is the critical density of the universe. Thus we get the density parameter 
Q = 0.6h x 10 n_12 . Setting this equal to unity (the QSSC does not have any limit 
on baryonic matter either from deuterium abundance or from CMBR anisotropy) 
we get for h = 0.6, a typical mass as 1.5 x 1O 13 M , suitable for a cluster. 
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Of course, as the above exercise shows, the results can be scaled up/down by 
rescaling the simulation parameters and thus are independent of the 'absolute' 
size of the box. A more detailed dynamical theory of the creation process will 
tell us how to relate the absolute size of clustering to the theoretical parameters. 




-2 -1.5 -1 -0.5 

log(r) 



Figure 6.8: The two point correlation function for the QSSC based model. Here, 
N = 100,000 and x = 0.3. As "time" goes on, the curve approaches more and more 
closely the slope of —1.8. The solid curve shows the result after 10 iterations. 

6.5 Results and Summary 

It must be stressed that these results are to be viewed as a preliminary report 
on a new scheme for generating structures in the Quasi Steady State Cosmology 
and quite a lot of follow up work has to be put into refining the model, so as 
to arrive at the values of the various parameters (chosen so far in an empirical 
way) from a deeper theoretical standpoint. Whatever the details of the creation 
process, the QSSC has repeated oscillations. We are trying to understand with 



the help of probability theory and stochastic processes, how clustering develops 
through such an iterative programme. 

The primary statistical indicator we have used in our analysis is the two point 
correlation function £(r). In order to further examine the statistical properties of 
the particle distribution, one must investigate the behaviour of higher moments 
and other quantities such as the "shape statistics". Work analyzing the higher 
moments, scaling relations , shape statistics etc using algorithms such as counts 
in cells, is in progress. 

However, the problem of formation of large scale structure being a complex 
one, it is desirable to keep the theoretical options open in the underlying cosmol- 
ogy. At the risk of stating the obvious, we should contrast the present approach 
from the standard approach to structure formation in the big bang cosmology. In 
the standard approach primordial fluctuations are postulated to begin with and 
their growth is studied under the effect of the gravitational field. Here, the main 
process which generates structures in the universe is the creation of matter around 
MCEs rather than gravitational instability. Our computer simulations show very 
clearly that the filament-cluster-void pattern observed in large scale structure 
can be generated simply from a creation algorithm. To what extent gravitational 
effects will further influence this picture remains to be seen. Although we have 
given preliminary ideas in subsection 6.2.1, a more detailed cosmogonic theory is 
also needed to tell us how coherent objects are ejected by mini-creation events. 

The success of the present toy model, however, holds out hope for a better 
understanding of structure formation via this alternative route. 



Chapter 7 
Conclusions 



I may not have gone where I intended to go, but I think I have ended up where I intended to be - Douglas 
Adams 

In the previous chapters we have established some new and important results 
in the field of nonlinear gravitational clustering and have shown that it is possible 
to gain insights into the role of gravity in clustering and consequent formation 
of large scale structure. The coherent theme presented in this work dealt with 
existence and models for the universal aspects of gravitational clustering, via 
various mechanisms. We also critically examined the role of gravity in structure 
formation theories by studying alternate models for structure formation, thus 
taking an open minded approach to the problem. 

We are led to the conclusion that there are universal phenomena in grav- 
itational clustering and we have analysed the reason behind their existence to 
some extent, analytically and semi analytically. In chapter 2 we established the 
existence of approximately invariant profiles for two point correlation function 
which evolved in a manner similar to linear growth, at all scales. It was also 
shown that this universal form is related to the fixed point phenomena in power 
transfer, which is already well known in literature. The connection between these 
forms "units of nonlinear universe" and the universal nonlinear scaling relation 
proposed by Hamilton and others is also clearly established. Numerical verifica- 
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tion of these results required high resolution N body simulations, which in turn 
led to analysis of two dimensional gravity in the next chapter. 

The next chapter tried to understand this result better by trying to under- 
stand the analytic framework underlying 2D simulations. In this work we have 
concluded that two dimensional structure formation simulations in an expanding 
background requires that we model it as parallel "infinite needles" with particles 
defined by the intersection of a plane orthogonal to these needles. We developed a 
(D+l) dimensional model of structure formation and have discovered that all the 
other approaches to 2D gravity violate either physical or consistency constraints. 

In the next chapter we have probed this connection further by addressing the 
exitence of such Nonlinear Scaling Relations in two dimensions. Two dimensional 
simulations provided us with much higher resolution so that this question could 
be addressed in detail. We have conclusively shown that although such a relation 
does exist in two dimensions as well, the behaviour at nonlinear end indicated 
that the clusters do not display the stable clustering behaviour as expected from 
3D simulations. Thus the slope at the nonlinear end for the scaling relation is 
not unity as would have been expected from a model based on standard stable 
clustering behaviour, contrary to some of the results claimed in literature. On the 
other hand, the results obtained supported the existence of a generalised 'stable 
clustering' that had been predicted earlier. 

We followed the trail of results and attempted to understand stable clustering 
and late time behaviour of clustering systems in the universe in more detail. 
This led to a detailed analysis of the behaviour of single cluster without making 
the usual additional assumptions about strict spherical symmetry and so on. A 
careful modelling of the asymmetries that are generated during the collapse phase 
of a density perturbation allowed us to derive a completely natural (as opposed to 



the ad hoc) procedure for stabilisation of a spherically collapsing system which was 
consistent with the usual reults. However our results differed quantitatively from 
the values obtained in the standard spherical collapse model. This allowed us to 
establish the connection between the shear and vorticity terms in the equations 
and the averaged pairwise velocity ratio function which is a measure of stable 
clustering. 

A critical examination of the role of gravity led us to investigate structure 
formation in Quasi Steady State Cosmology. The results of this investigation into 
a toy model without gravity (just expansion and creation of matter) indicate that 
it is possible to model some quantitative aspects of structure formation, such as 
the index of the two point correlation function by this process which is consis- 
tent with the accepted values. The visual images also reveal a clear cluster and 
void network which develops from a uniform distribution of particles. Statistical 
indicators such as minimal spanning trees which examine essentially all moments 
of the distribution reveal the growth of clustering in a manner comparable to 
growth of clustering in standard scenarios. Thus we have established that the 
observed two point correlation function with a slope of —1.8 can be obtained even 
in a model without gravitationally induced clustering. 

The analysis of dark matter gravitational clustering process has led to the 
following results established in this thesis. There exist generic aspects to gravita- 
tional clustering the most notable being the existence of fixed points in the way 
power is transferred from large scales to small scales. This transfer of power leads 
to scaling laws for statistical quantities such as two point correlation functions 
which are to a great degree independent of background cosmology and the matter 
power spectra. These scaling laws whose existence enable us to study late time 
behavior of the system easily have been vindicated in numerous simulations both 



in three and two dimensions. These two dimensional simulations which were in- 
tended to explore generic behavior of gravitational clustering led to a complete 
theoretical analysis of gravity in two dimensions which leads to greater coherence 
in the theoretical framework. 

It remains an open problem to explain and understand how the actual particle 
motion in the simulation volume is connected with power transfer and with the 
existence of fixed points in evolution of these quantities. Some work which ad- 
dresses the issue has indicated that the growth of power in small scales are caused 
by matter flows from large distances, essentially the central peaks gathering in 
the matter particles which flow along the regions of shocks or 'caustics' in the 
initial field. Thus the final configurations are connected to the initial Gaussian 
random field which are generated in inflationary scenarios for example. 

It can be observed that the final structures that are formed at late times 
tend to be spherical (when they are not tidally disrupted or have significant 
amounts of angular momentum). The existing spherical collapse models suffer 
from the ad hoc nature of the mechanism used to halt the spiraling collapse. A 
more rigorous analysis of spherical collapse model, including the amplification of 
asymmetries in the system in the collapsing phase leads to different thresholds 
for formation of collapsed structures with respect to the standard spherical model 
upon which theories such as the Press Schechter formalism is based. A study of 
this 'virialisation' process is important in modeling the formation of clusters of 
galaxies and in analyzing quantities such as mass functions and their evolution. 

These mass functions provide a quantitative measure of the dark matter wells 
in which the baryons can settle, cool and fragment and form the myriad structures 
that are observed. Connecting the observations with the semi empirical models of 
galaxy formation available currently as well as simulating the complete formation 



scenario from dark matter to clusters and galaxies is the fundamental program 
of structure formation theories in the standard model. 

This work is expected to lead us to a comprehensive model for baryonic 
structures. It is necessary to understand further the transfer of power in detail 
to explore generic features of gravity. This involves understanding second order 
effects which may be modelled by approximation schemes such as Zeldovich ap- 
proximation. It is possible to separate the particles into categories based on their 
net displacement and compute the contributions to the power spectrum from each 
category. This avenue must be explored to understand the details of connections 
between particle motion and power transfer. 

Subsequent to the formation of these clusters at high peaks at early times, 
the baryonic structures form a population of super massive stars, which serve 
as nuclei for galaxies at late times. A complete model for structure formation 
should be able to compute the initial star formation as well as the formation 
and propagation of radiation and ionization fronts in the neutral medium. This 
initial population of stars will lead to formation of super massive black holes and 
quasars and trigger of further star formation which may reionize the universe at 
a later epoch, leaving a signature on the cosmic microwave background. 

In parallel with this, an alternative cosmological structure formation scenario 
explored in this thesis is the quasi steady state cosmology. In studying this model 
and attempting to understand the observational data in the light of this model, 
we were led to a greater understanding of the standard model itself and a set of 
stringent tests that may be designed to distinguish between the many cosmologies. 
A large number of statistical measures, such as minimal spanning trees, fractal 
dimension, cluster mass functions, measures of shapes and orientations etc are 
used as discriminators between models which leads to a larger array of tests of the 



standard model per se which may cast further light on the formation processes. 
This model also brought to light an interesting question regarding power law 
profiles for two point correlation function which was generated even in a model 
devoid of gravitational clustering. 

Existence of other universal aspects of gravity, such as density and poten- 
tial profiles and their functional forms is another avenue which requires further 
exploration and is connected to the conclusions arrived at in this thesis. 

As the observational data drives us toward the era of precision cosmology, it 
will be possible to analyse these universal aspects in greater detail by comparing 
simulation models convolved with observational and instrument effects with maps 
of the sky which is expected to expose much of the weaknesses as well as reveal 
the strengths of the current paradigms of cosmology and structure formation. 
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